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SUMMARY 


This  report  sumnarizes  the  second  phase  of  research  in  nonlinear  signal 
enhancement  of  underwater  acoustic  measurements.  The  first  phase, 1  in  which 
the  nonlinear  process  was  developed  and  digitally  implemented,  generated  the 
need  for  a  mathematical  model  of  the  reflection  of  transient  acoustic  signals 
from  the  ocean  floor.  This  model  is  necessary  both  for  verifying  the  signal 
enhancement  process  by  predicting  the  ideal  response  of  a  typical  ocean  floor 
and  for  interpreting  measurements  that  have  been  enhanced.  Of  course,  many 
continuous  wave  models  are  available  for  prediction  purposes  but,  in  this  work, 
a  model  was  needed  that  could  predict  a  time-varying  field  and  relate  the  time- 
domain  features  to  specific  propagation  mechanisms. 

Several  classes  of  transient-response  models  have  been  developed  and  tested 
ranging  from  a  simple  ray  model  to  a  numerical  integration  of  the  complete  field 
integral.  The  simple  model  assumes  that  reflection  at  the  ocean  floor  can  be 
described  by  the  plane  wave  reflection  coefficient  and  the  remainder  of  the  path 
spreading  and  direction  can  be  determined  by  ray  tracing  (geometrical  acoustics). 
While  this  model  is  only  valid  for  high  frequency  or  deep  water,  it  is  flexible 
and  inexpensive  to  run.  Multiple  propagation  paths  are  easily  handled  as  are 
arbitrary  source  waveforms. 

At  the  other  extreme  is  the  integration  of  the  field  integral  and,  while 
this  procedure  can  be  expensive,  all  of  the  wave  effects  (interface  waves, 
diffraction)  are  properly  accounted  for.  For  low  frequency  (below  100  Hz), 
shallow  water  (several  wavelengths  or  less)  or  source/receiver  near  the  ocean 
floor,  these  wave  effects  can  be  critical.  Several  examples  of  application  of 
this  model  to  actual  measurements  are  given  in  which  the  various  wave  effects 
can  be  seen. 

Of  intermediate  complexity  and  cost  is  an  asymptotic  model  based  on 
approximate  solution  of  the  full  wave  integral.  This  solution  accounts  for 
multiple  paths  in  the  ocean  floor,  trapped  propagation  through  the  ocean  floor, 
and  field  displacement  on  reflection.  As  an  added  feature,  the  asymptotic 
solution  breaks  up  into  a  number  of  terms  each  of  which  corresponds  to  a 
physical  mode  of  propagation  in  the  bottom. 

While  these  models  were  designed  and  developed  primarily  as  an  aid  in 
the  interpretation  of  ocean  floor  seismic  measurements,  they  will  be  valuable 
in  any  circumstance  involving  propagation  of  transients. 


INTRODUCTION 


Modern  systems  for  acoustic  detection  and  localization  of  submarines  and 
surface  ships  must  function  in  a  wide  range  of  environments.  Detection  may 
be  desired  at  very  long  ranges  if  large  ocean  areas  must  be  protected  or  at 
shorter  ranges  if  the  guard  areas  are  straits  or  passages.  Again,  depending  on 
the  area,  detection  and  localization  may  have  to  be  done  in  deep  or  shallow 
water,  over  rough  or  smooth  regions  of  the  ocean  floor  and  in  many  different 
conditions  of  water  temperature  structure.  Depending  on  the  specific  target 
and  tactics,  the  detection  or  localization  system  may  operate  across  a 
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wide  spread  of  frequencies  and  the  sensor  itself  (the  hydrophone  or  hydrophone 
array)  may  be  anywhere  in  the  water  column  from  near  the  surface  to  on  the 
bottom. 

In  some  of  these  operating  environments,  the  ocean  floor  plays  only  a 
minor  role  in  the  transmission  of  acoustic  energy  from  source  to  sensor  but 
there  are  three  cases  in  which  the  ocean  floor  can  greatly  influence  propagation: 

a.  low  acoustic  frequency  -  below  about  100  Hz  in  typical  ocean  situations, 

b.  shallow  water  -  acoustic  wavelength  is  of  the  order  of  water  depth,  and 

c.  near-bottom  sensors  -  sensor  is  within  an  order  of  the  wavelength  of 
the  bottom. 

For  sensor  environments  other  than  these  three,  the  ocean  floor  can  often  be 
adequately  treated  as  a  planar,  lossy  reflector.  To  do  so  in  the  above- 
mentioned  situations  will,  however,  introduce  serious  errors  into  the  predicted 
behavior  of  the  bottom. 

Since  accurate  prediction  of  acoustic  propagation  is  crucial  not  only  to 
design  and  evaluation  of  new  sensors  but  also  to  design  of  tactics,  knowledge 
of  the  acoustical  behavior  of  the  ocean  floor  is  a  vital  part  of  sonar  sensor 
research  and  application.  Unfortunately,  the  ocean  floor  has  still  not  been 
well  characterized  at  low  frequency.  Thousands  of  measurements  of  ocean 
bottom  properties  have  been  made  by  many  laboratories  but  most  of  these 
measurements  are  either  not  applicable  to  system  operation  below  100  Hz  or 
have  not  been  analyzed  for  this  information  in  particular.  The  NAVAIRDEVCEN 
itself  has  hundreds  of  seismic  measurements  from  sites  around  the  world2  and 
these  are  only  recently  being  enhanced*  sufficiently  to  be  useful  for  detailed 
bottom  property  estimation. 

Because  the  ocean  floor  is  often  responsible  for  several  distinct  paths 
of  propagation,  short-duration  pulses  are  ideal  for  making  the  necessary 
property  measurements.  A  sharply  pulsed  source  results  in  separate  received 
pulses  corresponding  to  each  path  from  source  to  receiver.  Transient  signals 
are,  however,  considerably  more  difficult  to  model  theoretically  than  continuous 
wave  (CW)  signals  and  this  theoretical  modeling  is  necessary  in  order  to 
interpret  measurements  physically. 

Assumption  of  continuous,  sinusoidal  excitation  causes  the  wave  equation 
to  degenerate  into  a  simpler  differential  equation  that  can  be  solved  by  (for 
example)  normal  modes  in  a  straightforward  manner.  If,  on  the  other  hand,  the 
excitation  is  transient  then  either  the  wave  equation  must  be  solved  directly 
or  the  CW  solution  must  be  integrated  over  frequency.  In  this  report,  we  will 
not  consider  direct  solution  of  the  wave  equation  but  rather  we  will  develop 
several  solution  techniques  based  on  integration  of  the  CW  solution.  While 
CW  solutions  are  common  and  highly  developed,  transient  or  time-dependent 
solutions  are  not,  so  the  models  developed  as  a  part  of  this  work  will  be 
generally  useful  as  an  adjunct  to  acoustic  studies  of  the  ocean  floor. 
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THE  PHYSICAL  MODEL 


Before  we  actually  construct  a  model  for  the  acoustical  behavior  of  the 
ocean  floor,  let  us  consider  the  relevant  acoustic  properties.  The  most 
important  property,  sound  speed,  is  the  speed  at  which  a  small  disturbance 
travels  through  the  medium.  Changes  in  sound  speed  within  the  medium  cause 
refraction  or  bending  of  the  sound  waves  (primary  effect)  and  diffraction  of 
the  field  (secondary  effect)  into  regions  not  reached  by  geometrical  rays. 
Discontinuities  in  sound  speed  occurring  at  boundaries  between  media  of 
different  types  cause  reflection  (primary  effect)  and  formation  of  interface 
waves  (secondary  effect). 

In  fluids,  acoustic  waves  are  compressional  waves  which  are  characterized 
by  a  compressional  sound  speed  that  is  a  function  of  both  the  density  and  the 
compressibility  of  the  fluid.  Thus,  we  will  also  need  to  know  the  density 
(or  the  compressibility)  of  the  fluid  in  order  to  describe  its  behavior  for 
no-loss  propagation. 

Elastic  solids,  on  the  other  hand,  support  two  distinct  types  of  waves: 
compressional  and  shear  waves.  We  then  have  another  sound  speed  based  on  the 
density  and  resistance  to  shear  or  twisting  deformation  of  the  solid  and 
another  set  of  waves  that  can  be  reflected  or  refracted  and  can  form  interface 
waves.  We  will  also  consider  dissipation  or  attenuation  of  shear  and  com¬ 
pressional  waves  in  order  to  represent  the  real  ocean  floor  materials  more 
accurately. 

Certain  simplifications  will  need  to  be  made  in  our  representation  of 
these  properties  so  that  the  resulting  model  is  not  prohibitively  complicated. 

For  example,  we  will  not  consider  boundary  roughness  or  slope.  Since  our 
principal  interest  is  in  low  frequency  propagation,  this  will  not  be  too 
restrictive  an  assumption.  In  addition,  we  will  restrict  variations  in  sound 
speed  to  one  dimension  (depth)  and  further  restrict  the  sound  speed  to  be 
either  piece-wise  constant  or  of  an  analytical  form  that  allows  solution  in 
terms  of  special  functions. 

One  final  assumption  concerns  the  manner  in  which  attenuation  is  modeled. 
Particularly  in  fluid-saturated  sediments,  the  mechanisms  for  energy  dissipation 
can  be  complicated.  The  usual  viscous  resistance  of  the  fluid  is  combined  with 
the  friction  between  sediment  grains  and,  possibly,  the  interaction  of  entrained 
gases.  Data  are  sparse  on  these  processes  and,  even  if  measurements  were 
available,  incorporation  of  such  a  complicated  process  into  a  computer  program 
would  be  a  major  undertaking.  Consequently,  we  will  approximate  any  attenuation 
process  as  a  viscous  damping  process  in  which  the  dissipation  per  cycle  is 
constant.  This  is  cnly  valid  over  limited  ranges  of  frequency. 

Figure  1  illustrates  the  physical  model  that  will  be  used  for  most  of  this 
research.  This  model  consists  of  a  solid  basement  half-space  over  which  lies  a 
fluid  layer  with  variable  properties.  This  fluid  layer  can  either  have  constant 
sound  speed  in  which  case  the  field  variables  will  vary  as  linear  combinations 
of  sines  and  cosines  or  1/c^  linear  in  depth  in  which  case  the  solutions  will 
be  linear  combinations  of  Airy  functions.  Since  many  models  are  available  for 
acoustic  propagation  in  water,  we  will  not  consider  this  part  of  the  problem; 
hence,  we  will  use  a  homogeneous  half-space  of  fluid  above  our  ocean  floor  model. 
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GEOPHYSICAL  MODEL 


•>> 
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FLUID 


Figure  1. 


Generalized  physical  model  of  ocean  floor. 
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While  this  is  a  very  simple  model,  it  does  allow  most  of  the  important 
types  of  acoustic  propagation  to  be  represented.  Furthermore,  addition  of 
layers  is  straightforward  and  would  permit  more  accurate  modeling  of  sound  speed 
profiles  in  sediments. 

In  the  following  two  sections,  we  will  discuss  the  two  major  theoretical 
approaches  toward  calculation  of  the  reflected  field  from  this  bottom  model. 

The  simpler  procedure  involves  assuming  that  the  incident  energy  is  in  the  form 
of  plane  waves.  The  plane  wave  reflection  coefficient  is  relatively  easy  to. 
compute,  and  a  fairly  complicated  field  can  be  approximated  by  adding  different 
plane  wave  solutions  together.  Unfortunately,  a  number  of  interesting  and 
sometimes  significant  effects  are  not  treated.  Interface  waves,  for  example, 
are  not  excited  according  to  plane  wave  theory. 

The  second  type  of  solution  procedure  takes  into  account  the  curvature  of 
the  waves  emanating  from  a  point  source.  This  method  is  quite  time-consuming 
but  it  does  provide  an  exact  solution  and  will  be  useful  in  establishing  the 
extent  of  the  effects  omitted  by  the  plane  wave  solution. 


PLANE  WAVE  SOLUTION 


Although  it  exists  in  three-dimensional  space,  a  plane  wave  is  essentially 
one-dimensional.  Its  wavefront  is  a  plane  and  the  medium  is  displaced  by  its 
passage  only  in  the  direction  perpendicular  to  this  plane.  A  true  plane  wave 
can  only  remain  plane  in  a  medium  of  constant  sound  speed;  however,  we  will 
consider  waves  that  are  approximately  plane  (or,  locally  plane)  so  that  we 
can  apply  this  theory  to  our  inhomogeneous  ocean  floor  model. 


True  plane  waves  are  solutions  of  the  three-dimensional  wave  equation 
for  constant  sound  speed, 


32»  +  d^±  +  _  1  32i> 

ax2  ay2  sz2  c 7  It? 


(1) 


where 


c  =  sound  speed 

p  =  some  scalar  property  of  the  field  (pressure,  velocity  potential). 

If  we  assume  harmonic  functions  of  the  form  e”1;jt  and  separate  the  field 
quantity  <p  into  a  range-dependent  factor  and  a  depth-dependent  factor,  we  can 
write  three  ordinary  differential  equations:  two  for  the  range  function  and 
one  for  the  depth  function.  The  depth  function  equation  is, 

+  Ck2(z)-Y2]  U  =  0  (2) 

dz^ 
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where 

k  *  u/c(z) 

U  »  depth  function  of  field  property 
y  «  horizontal  wave  number 

If  the  sound  speed  is  constant,  the  solutions  of  equation  (2)  are, 
U  -  ei'62 

where 


B  *  vertical  wave  number 
= Vk2-y2'  so  that  k2  =  y2+e2. 

In  conjunction  with  the  harmonic  factor  e"1ut,  these  solutions  represent  waves 
with  a  vertical  component  of  motion:  the  positive  sign  corresponds  to  movement 
in  the  +z  direction  and  the  negative  sign  to  movement  in  the  -z  direction.  This 
is  true  in  general  for  the  dependence.  If  the  coefficient  of  a  wave 

number  component  is  positive  then  the  wavefront  is  moving  in  the  positive 
direction  along  the  corresponding  axis.  Keep  this  feature  in  mind  when  we 
discuss  the  equations  for  reflection.  This  distinction  between  positive  and 
negative-going  waves  will  allow  separation  of  the  incident  and  reflected  fields. 

The  complete  solution  to  equation  (1)  is. 


*  =  *oe 


i[k-R-wt] 


where 


k  =  vector  wave  number,  y  i  +  y  j  +  Sk 

a  y 

R  =  position  vector,  xi  +  yj  +  zk 
uj  =  angular  frequency,  2nf 
$o  =  amplitude 


Yx2  +  YyZ  *  S2  =  k2  =  (to/c)2 . 

The  vector  wave  number  gives  the  direction  of  the  normal  to  the  wavefront  and 
jo  gives  the  frequency  of  oscillation  for  the  single  frequency  wave  or  equation  (3). 

Since  the  properties  of  a  plane  wave  are  identical  anywhere  on  a  plane 
perpendicular  to  the  wave  number  vector  (i.e.,  on  a  wavefront),  boundary 
conditions  are  very  simple  to  write.  For  a  fluid-fluid  boundary,  we  must  maintain 
continuity  of  displacement  (or  velocity)  normal  to  the  boundary  and  continuity 
of  pressure  across  the  boundary.  At  a  fluid-solid  boundary,  the  condition  of 
zero  resistance  to  transverse  displacement  of  the  solid  is  also  needed. 
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Only  the  lower  half-space  of  our  physical  model  is  elastic  so  we  will  not 
consider  propagation  in  elastic  media;  just  reflection  from  them.  The  plane 
wave  reflection  coefficient  is  the  ratio  of  reflected  field  amplitude  to 
incident  field  amplitude.  For  a  fluid  over  an  elastic  half-spaced,  the  plane 
wave  reflection  coefficient  is, 3 


V  = 


G^-G 

Gj+G 


(4) 


where 


G  = 


G 


1 


B 


1,2 


B 


s 


B  = 


c 


1 


C2 


c 


V 


f  = 


“1 


"2 


p1^2//p2 

B^Uf-l)2  +  4f2B2Bs) 

=  V(cv/Ci>2)  ■  1 

shear  sound  speed  of  solid 
compressional  sound  speed  of  fluid 

compressional  sound  speed  of  solid 

vertex  velocity  of  wave  (y=io/cv) 

62/C2v 

density  of  fluid 
densi ty  of  sol  id 


Because  we  will  be  placing  a  layer  over  the  solid,  we  will  need  to  use 
this  reflection  coefficient  as  a  boundary  condition. 4  First,  we  write  the 
pressure  in  a  homogeneous  fluid  overlying  the  solid  half-space  (the  z-axis  being 
perpendicular  to  the  interface). 


P 


+  Ve 


where 


PQ  =  pressure  amplitude  of  incident  wave 
i  =  vertical  wave  number 
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and  the  normal  derivative  of  the  pressure  (which  is  proportional  to  particle 
velocity) , 

H  -  tePo  [-*'UZ  *  ve'8z] 

The  ratio  of  these  quantities  is  proportional  to  acoustic  impedance  and  can 
be  used  as  a  boundary  condition.  At  the  interface,  z=0;  therefore  this 
impedance  is. 


-  p  -  i ( 1+V) 

o  ■  3p/3Z  Z=Q  e(l-V) 


(5) 


This  relation  will  be  identical  for  any  quantity  that  is  a  solution  of  the 
wave  equation. 

At  this  point,  we  must  also  consider  the  addition  of  attenuation  to  this 
model.  We  will  assume  viscous  attenuation  which  is  a  constant  loss  per  cycle 
equivalent  to  the  addition  of  imaginary  part  (linear  in  frequency)  to  the 
wave  number, 

k'=k+iwa  (6) 

so  that  the  plane  wave  function  becomes, 

ik'R  -waR  ikR 
e  =  e  e 

and  the  loss  A  (in  decibels)  over  some  distance  A r  is, 

A  -  -20  log  [e-"'r] 


=  8.69  waAr 

This  loss  is  generally  given  in  decibels  per  distance  Ar  at  some  frequency  f, 
therefore, 


-  A 

a  "  54.6  fAr 

In  practice,  the  attenuation  is  inserted  by  modifying  the  sound  speed,  not 
the  wave  number.  Equation  (6)  can  be  written, 

k‘  =  oj/C1  =  w[l  +  iac)/c 


therefore. 


c1  =  c/(l  +  iac)  (8) 

This  substitution  will  allow  the  ad  hoc  addition  of  small  amounts  of  attenuation 
to  relations  derived  for  conservative  systems  like  the  reflection  coefficient 
equation  (4). 
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Next,  let  us  consider  reflection  of  plane  waves  from  a  fluid  layer  over 
a  solid  basement.  Appendix  A  presents  a  solution  for  the  general  case  of  a 
multilayered  reflector,  so  we  can  just  simplify  these  results.  For  a  constant 
sound  speed  layer,  the  model  is  shown  in  Figure  2.  In  the  fluid  half-space, 
the  field  can  be  written  as, 

-i6  z  is  z 

*o  3  e  +  vte 


where 


Vt  =  plane  wave  reflection  coefficient  of  layer  plus  basement 


=  velocity  potential. 


Also,  the  field  in  the  layer  is, 

<|>j  =  a1  cosSjZ  +  b^  sin^z. 

Using  equation  (A-4) ,  we  can  evaluate  the  ratio  of  coefficients  of  <#>- , 


sine^  +  4  cosb^ 
1  ”  cosB^  -  5  sinBjh 


r,  = 


(9) 


where  t,  is  given  by  equations  (4)  and  (5).  The  overall  reflection  coefficient 
is  then  given  by  equations  (A-2)  and  (A-3) , 


_  !  8l  +  ir^/^ 

t  rQ  -01  +  ir1B0P1/p0 


(10) 


We  will  also  consider  a  layer  with  variable  sound  speed:  in  particular, 

l/c2(z)  =  1/c2  -  gxz  (11) 

This  model  is  shown  in  figure  3  and  the  solutions  associated  with  this  sound 
speed  function  are  discussed  in  appendix  B.  The  resulting  field  expression 
is  a  linear  combination  of  these  solutions  (called  Airy  functions). 


=  ajA-Uj)  +  b1Bi(Z1) 

where 

A.. ,  B..  =  Airy  functions 
and,  from  equation  (B-5), 

2 


(12) 


l1  =  qz 


■?(H 


(13) 
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FLUID  p 

cO  ° 

FLUID 

SOLID  pB  b 

S' 

CB 

Figure  2.  Single,  constant  sound  speed  layer  over  solid  basement. 
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where 


q  =  (9j«  ) 


Therefore , 


5T-  ’  5T  fi  ‘  <z»-  |14) 

As  in  the  case  of  the  constant  sound  speed  layer,  we  solve  for  the  overall 
reflection  coefficient  by  means  of  equations  (A-4),  (A-3),  and  (A-2)  of 
appendix  A.  These  yield,  respectively, 

VZ1>-«VZ1>  ^ 

-P0[r1Aj  +  Bl]  ♦  <a0Pt  IrtA,  »  B,1  (16) 

r°  »„lrlAi  +  V  *  'Vl  |rlAi  +  B(J 


Vt  -  l/r0. 

We  now  are  able  to  compute  the  reflection  coefficient  for  either  a  constant 
sound  speed  layer  over  a  solid  (equations  (5),  (9),  and  (10))  or  a  gradient 
layer  over  a  solid  (equations  (5),  (15),  and  (16)).  Each  of  these  equation 
sets  provides  a  value  for  for  a  specific  frequency  and  a  specific  incident 
angle  of  the  wavefront.  This  incident  angle  is  actually  determined  by  the 
components  of  the  wave  number  vector  so  that  we  will  use  y»  the  horizontal 
wave  number  as  the  independent  variable.  The  relationships  between  y»  6>  and 
the  grazing  angle  9  (angle  to  the  horizontal  plane)  of  the  wavefront  normal  are, 

y  =  k  cose 

6  =  k  sine  (17) 

2  2  2 

Although  these  single-layer  models  are  very  simple,  they  do  include  most  of 
the  important  types  of  propagation.  The  gradient  layer  model  allows  reflection 
from  the  top  of  the  layer,  refraction  within  the  layer,  reflection  from  the 
solid  basement,  and  interface  propagation  along  both  the  fluid-fluid  and  the 
fluid-solid  interfaces. 
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As  we  have  mentioned,  one  calculation  of  Vt  is  needed  for  each  frequency 
and  each  grazing  angle.  In  order  to  compute  the  response  of  this  model  to 
an  impulse,  we  must  first  compute  the  reflection  coefficient  at  many  frequencies 
throughout  the  band  of  interest.  The  inverse  Fourier  transform  then  yields 
the  impulse  response  in  the  time  domain.  For  a  specific  incidence  angle 
(which  we  would  normally  compute  for  a  particular  path  by  ray  tracing  through 
the  ocean  medium),  the  time-domain  response  is  as  follows, 


p(t)  =  ^  j  Vt(w)e"1ujtdoj  (18) 


where 


p  =  acoustic  pressure. 

(The  -iwt  exponent  reflects  our  choice  of  e"iut  as  the  harmonic  time  dependence 
of  the  acoustic  disturbance.  This  is  contrary  to  the  convention  in  signal 
processing.)  While  we  have  derived  Vt  in  terms  of  the  velocity  potential  <f>, 
we  can  consider  it  to  be  the  pressure  reflection  coefficient  also  since  both 
p  and  <j>  identically  satisfy  the  wave  equation.  Equation  (18)  gives  the  pressure 
time-history  for  the  reflection  of  a  unit  pressure  impulse.  This  result  is 
identical  to  the  velocity  potential  time-history  for  an  incident  unit  impulse 
of  velocity  potential.  A  source  of  unit  pressure  impulse  is  considerably 
easier  to  understand  than  a  unit  impulse  of  velocity  potential  so  we  will 
continue  our  discussion  in  terms  of  pressure  only. 

Once  we  have  found  the  impulse  response,  we  can  add  a  source  waveform  in 
one  of  two  ways.  In  the  time  domain,  the  source  time  function  can  be  convolved 
with  the  impulse  response  to  produce  the  time  history  for  reflection  of  that 
source  waveform, 

r(t)  =  p(t)  *  s(t)  (19) 


where 


p(t)  =  impulse  response  of  reflector  by  equation  (18) 
r(t)  =  reflected  signal 
s(t)  =  source  waveform. 

Alternatively,  we  can  make  use  of  the  fact  that  a  convolution  in  the  time  domain 
is  a  multiplication  in  the  frequency  domain  and  multiply  the  reflection 
coefficient  Vt  by  the  source  spectrum  prior  to  transforming  so  that. 


Vt(w)  S(w)e'1ujtdw 


(20) 
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where 


$(w)  *  Fourier  transform  of  s(t) 
s(t)eiwtdt. 


/ 


In  short,  the  procedure  for  computing  the  time-domain  response  for  a 
particular  propagation  path  given  a  nonsinusoidal  source  is  as  follows: 
compute  the  water  path  from  source  to  ocean  floor  to  receiver  by  ray  tracing. 
Using  the  grazing  angle  (relative  to  the  ocean  floor  at  the  point  of  reflection) 
of  this  ray,  calculate  the  horizontal  wave  number  and  then  compute  Vt  at  many 
frequencies  across  the  band  of  interest.  Adjust  this  result  for  the  spreading, 
time  delay,  and  phase  shifting  incurred  along  the  ray  path  in  the  water. 

Multiply  this  function  by  the  spectrum  of  the  source  waveshape  and  take  the 
inverse  Fourier  transform.  The  result  is  then  the  received  signal  for  the 
selected  propagation  path.  Plane  wave  reflection  theory  is  used  to  account 
for  the  effects  of  the  reflection  model  (the  layer  and  solid  basement,  in  this 
case)  while  ray  theory  is  used  to  predict  the  water-borne  effects.  We  can 
repeat  this  procedure  for  as  many  paths  as  are  significant  in  a  particular 
environment  and  combine  the  results,  as  shown  in  figure  4,  to  produce  a  complete 
prediction  for  the  received  signal. 

This  technique  for  building  a  signal  from  its  component  parts  is  a  flexible 
and  computationally  fast  method  for  predicting  the  propagation  of  transient 
acoustic  energy.  It  is,  however,  limited  by  the  restrictions  of  both  ray 
theory  and  plane  wave  reflection  theory.  In  instances  in  which  the  source 
or  receiver  is  near  the  ocean  floor,  or  in  shallow  water,  or  at  low  frequency, 
this  approach  may  be  inadequate.  In  these  circumstances,  interface  waves, 
both  compressional  and  shear,  may  contribute  substantially  to  the  field. 

There  may  also  be  considerable  diffraction  at  low  frequency  which  cannot  be 
predicted  by  ray  tracing  methods.  For  these  reasons,  we  will  consider  in  the 
next  section  a  more  sophisticated  and,  consequently,  a  more  expensive  model 
that  represents  the  wave  field  itself  without  the  approximations  implicit  in 
ray  tracing  and  plane  wave  reflection. 


SPHERICAL  WAVE  SOLUTION 


A  small  source  (small  in  relation  to  the  dimensions  of  the  medium  and 
boundaries)  generates  a  spherical  wavefront  in  a  homogeneous  medium  and  such 
a  wave  acts  differently  at  boundary  reflection  than  a  plane  wave.  Usually, 
the  assumption  is  made  that  the  source  is  far  enough  removed  from  any  boundaries 
so  that  the  wavefront  may  be  considered  plane.  As  we  will  see,  this  assumption 
cannot  be  justified  at  low  frequency  or  when  the  source  is  near  the  ocean  bottom. 
We  will  develop  a  quantitative  measure  by  which  we  can  determine  the  validity 
of  the  plane  wave  assumption  in  any  practical  case  but  for  now  it  is  sufficient 
to  recognize  that  in  some  instances  wavefront  curvature  cannot  be  neglected. 
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Figure  4.  Ray/plane-wave  composite  model  for  transient  propagation. 
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In  a  homogeneous  fluid,  the  free-field  expression  for  a  harmonic  point 
source  is. 


fii (kR-wt) 
♦  =  - R 

where 


(21) 


k  =  wave  number,  w/c 
R  =  distance  from  source 
<p  =  velocity  potential. 


This  is  the  simplest  expression  of  a  spherical  wave  and  it  is  true  not  only 
for  velocity  potential  but  for  any  field  quantity  (pressure,  particle  velocity, 
displacement  potential)  that  satisfies  the  Helmholtz  equation  for  acoustic 
vibrations.  Also,  the  expression  is  true  in  the  immediate  vicinity  of  a 
point  source  in  an  inhomogeneous  medium. 

Unfortunately,  equation  (21)  is  spherically  symmetric  while  the  boundaries 
are  typically  parallel  planes.  Therefore,  either  cylindrical  or  Cartesian 
coordinates  would  be  more  suitable  for  ocean  propagation  problems.  Just  as 
the  spherical  wave,  given  by  equation  (21),  is  the  elementary  solution  of  the 
vibration  equation  in  spherical  coordinates,  the  plane  wave,  equation  (3),  is 
the  elementary  solution  in  Cartesian  coordinates  and  the  cylindrical  wave, 

*  *  J0(yr)ei(-ez_wt)  (22) 


is  the  elementary  solution  in  cylindrical  coordinates. 

Since  an  integral  superposition  over  wave  number  of  the  elementary  plane 
wave  solutions  is  a  Fourier  transform  integral  and  the  result  must  be  identical 
to  the  free-field  expression,  equation  (21),  we  can  invert  the  transform  to 
find  the  weighting  function. 5  The  complete  superposition  integral  can  then  be 
written. 


1 (kR-wt) 
- 


ik-R 
e _ 

6 


dy  dy 
x  y 


(23) 


where 


y  components  of  horizontal  wave  number  y. 


Notice  that  the 
and  a  weighting 


i  k 

inteqral  comprises  the  elementary  plane  wave  expression  e 
tunction  i /  (2-r 3 ) . 
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Furthermore,  an  integral  superposition  of  the  elementary  cylindrical  waves 
is  a  Hankel  transform  and  can  be  inverted  similarly.  The  result  is  identical 
to  a  conversion  of  equation  (23)  to  cylindrical  coordinates. 


i (kR-wt) 


i  e-'"'  j 


J„(vr)ei162 


ydy 


(24) 


Figure  5  illustrates  the  general  reflection  geometry  that  we  will  consider. 
In  particular,  a  point  source  and  receiver  are  located  above  an  inhomogeneous 
half-space  that  is  characterized  by  a  plane  wave  reflection  coefficient. 

Since  equations  (23)  and  (24)  were  derived  for  the  source  located  at  z=0,  we 
will  have  to  modify  the  expressions  accordingly. 

Actually,  two  field  expressions  are  required:  one  for  the  propagation 
directly  from  the  source  to  the  receiver  (r  units  horizontally  and  d-z  units 
vertically);  and  one  for  the  propagation  from  the  source  down  to  the  half-space, 
modification  by  the  reflection  coefficient,  and  up  to  the  receiver  (r  units 
horizontally  and  a  total  of  d+z  units  vertically).  Only  this  second  field  is 
of  interest  here.  In  terms  of  plane  waves,  this  field  is, 


i(yx*  +  yy y)  is(z+d) 
e _ *  e _ 

6 


v(y)  dyxdyy 


or,  in  cylindrical  coordinates, 

i6(z+d) 

-  V(y)ydy.  (25) 

o 

This  integral  can  be  interpreted  as  follows:  first,  the  point  source  is 
represented  as  an  infinite  sum  of  plane  waves  each  travelling  in  the  direction 
given  by  the  horizontal  wave  number  y.  The  horizontal  wave  number  is  related 
to  the  constant  wave  number  by  the  angle  of  the  direction  of  propagation  with 
the  horizontal, 


. .  J  ^ 


J„(yr)e 


y  =  kcoS1. 

Hence,  as  ,  increases,  the  plane  wave  travels  more  nearly  horizontal.  Eventually, 
,  is  greater  than  k  and  this  angle  becomes  imaginary.  These  waves  are  evanescent 
waves.  They  travel  horizontally  but  the  field  amplitude  decays  exponentially 
in  the  vertical  direction.  This  can  be  seen  by  substituting  an  imaginary  &  into 
the  elementary  plane  wave  expression,  equation  (3).  The  vertical  wave  number, 
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becomes  imaginary  when  y  is  greater  than  k.  The  integral  of  equation  (25)  thus 
states  that  a  point  source  cannot  be  completely  described  by  superposition  of 
real  plane  wave  but  these  evanescent  waves  must  be  included  also. 

The  second  point  to  notice  in  equation  (25)  is  that  the  elementary  waves 
inside  the  integral  are  each  modified  by  the  appropriate  value  of  the  plane 
wave  reflection  coefficient.  In  this  way,  the  influence  of  the  inhomogeneous 
reflector  is  accounted  for  completely.  Finally,  the  exponent  involving  the 
vertical  wave  number  is  positive  which  correctly  describes  the  upward  direction 
of  travel  of  the  reflected  field. 

Equation  (25)  then  describes  the  field  at  the  receiver  po'nt  including 
all  of  the  effects  of  the  reflecting  medium.  Unlike  the  plane  wave  model, 
this  expression  explicitly  involves  the  source  and  receiver  locations  so  that 
the  solutions  will  not  be  as  simple  to  extend  to  other  geometries  as  they  were 
with  the  plane  wave  model.  We  do,  however,  have  a  means  for  describing  all 
of  the  types  of  propagation  that  result  from  the  reflecting  medium.  Before 
we  solve  this  equation  numerically,  we  will  examine  some  of  the  properties 
of  this  equation  through  an  approximate  analytical  technique. 

By  using  the  identities,6 

J0(yr)  =  i  [H0{l)(Yr)  +  H0(2)(Yr)j 

H0(1)(yr)  =  -Ho(2)(-Yr) 

where 

H^,  HQ'2^  =  Hankel  functions, 
we  can  rewrite  equation  (25)  in  the  following  form, 


•  -iwt  f  ie(z+d) 

i  j -  J  H0  1  (vr)  - — -  V( y )  ydy  (26) 

— oo 

as  long  as  V(y)  =  V(-y).  This,  of  course,  is  always  true  for  physically 
realizable  problems  since  it  does  not  matter  whether  the  plane  wave  is  coming 
from  the  right  or  the  left  -  only  the  magnitude  of  the  angle  is  important. 

If  <r  is  large  compared  to  1,  we  can  further  reduce  the  integral  by  using 
the  asymptotic  form  of  H0(l), 
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Equation  (26)  becomes. 


.  -i(uit  +  ir/4) 

♦  =  - 

yjzTP 


This  is  the  form  of  the  integral  that  we  will  consider  in  the  following 
approximate  techniques.  Although  we  will  not  consider  this  aspect  until  the 
next  section,  notice  that  the  integral  itself  has  the  form  of  a  Fourier  transform. 


I  =  J  F(y)eiYr  dY  (28) 


where 

F(y)  =  ltd  A  ei6(z  +  d) 

As  a  result,  the  solution  can  be  computed  by  means  of  the  fast  Fourier  transform 
(FFT)  algorithm. 

We  can  learn  quite  a  bit  about  the  physics  of  reflection  by  considering 
the  integral  in  equation  (27)  as  a  complex  contour  integral.  In  general, 
the  reflection  coefficient  V  will  have  poles  for  some  complex  (and  possibly, 
some  real)  values  of  y  and  a  nontrivial  branch  cut  emanating  from  the  value 
of  v  at  which  3=0.  As  we  will  see,  each  of  these  singularities  is  related  to 
a  specific  type  of  propagation  and  it  is  possible  to  isolate  these  mechanisms 
in  the  solution  process. 

Given  ttiat  equation  (27)  is  a  complex  contour  integral,  there  are  several 
ways  to  find  a  solution.  First,  the  integral  can  be  evaluated  (numerically) 
along  the  original  contour  which  is  the  real  axis  of  the  y-plane.  Second,  the 
contour  can  be  deformed  in  such  a  way  that  only  pole  residues  and  branch  line 
integrals  remain  and  the  branch  line  integrals  can  be  evaluated  numerically. 
Also,  the  contour  can  be  deformed  so  that  it  passes  through  the  saddle  points 
of  the  integrand.  The  contribution  to  the  integral's  value  becomes  concen¬ 
trated  at  these  saddle  points  and  the  integral  can  be  evaluated  approximately 
in  the  neighborhood  of  each  saddle  and  the  remainder  of  the  path  can  be 
neglected.  Finally,  various  combinations  of  these  methods  may  be  used  as  long 
as  the  rules  for  contour  deformation  are  observed. 

The  direct  integration  method  and  the  residue/branch-1 ine-integral  solution 
are  standard  methods  so  we  will  only  consider  them  briefly  in  subsequent 
discussions.  The  steepest  descent  method  or  saddle  point  method  is  based  on 
locating  the  regions  in  the  .-plane  in  which  the  integral  becomes  concentrated 
and  then  determining  the  integral's  value  in  these  regions. 3>7  since  the 
integral  is  concentrated  in  a  small  area,  the  value  can  be  found  by  series 
approximation.  In  addition,  the  regions  of  concentration  (the  saddle  points) 
correspond  to  physical  propagation  mechanisms. 


/ 


,i  [yr  +  e(z  +  d)] 

3 


V(y)  V Y  dy 


(27) 
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To  locate  the  saddle  points,  we  first  separate  the  integrand  into  a  slowly 
varying  magnitude  function  and  a  complex  exponential.  In  our  case,  the 
reflection  coefficient  phase  can  vary  quite  rapidly  so  we  will  have  to  consider 
this  phase  as  part  of  the  complex  exponential.  Thus,  we  can  write  the  integral 
part  of  equation  (27)  as. 


I 


ei  hr  +  s(z  +  d)  +  x)d 


(29) 


where 

V(y)  =  Vm(Y)e1?(y) 

and  we  will  assume  that  the  factor  Vm/7/S  is  slowly  varying.  This  is  not,  in 
fact,  true  in  some  instances  (particularly  near  8=0)  but  as  long  as  we  are 
aware  of  the  assumption  we  will  be  able  to  compensate  for  any  problems  it 
introduces. 

Since  the  exponent  of  the  exponential  in  equation  (29)  is,  in  general, 
complex,  it  controls  the  rate  of  change  not  only  of  the  phase  of  the  integrand 
but  also  of  the  magnitude  (as  long  as  the  rest  of  the  integrand  is  slowly 
varying).  We  would  expect  that  if  the  phase  is  changing  rapidly,  the  individual 
elements  of  the  integration  will  add  out  of  phase  and,  therefore,  contribute 
little  to  the  integral's  value.  On  the  other  hand,  when  the  phase  of  the 
integrand  is  stationary,  the  elements  will  add  constructively.  At  the  saddle 
point,  both  the  magnitude  and  phase  are  stationary  and  we  can  find  a  path 
through  this  point  along  which  the  magnitude  decay  is  most  rapid.  This  is 
the  path  of  steepest  descent  and,  by  virtue  of  the  behavior  of  analytic  functions, 
the  phase  along  this  path  is  constant.  Physically,  this  means  that  the  integral 
is  actually  made  of  discrete  groups  of  wavefront  arrivals.  Each  group  is  given 
by  a  small  amount  of  waves  of  slightly  different  wave  number  that  are  in-phase 
at  the  receiver  point.  The  saddle  points  locate  each  of  these  groups. 

In  order  to  locate  the  saddle  points  and  descent  paths,  rewrite  equation 
(29)  in  the  following  form. 


I 


/ 


G  (  y  )  e 


f(-.) 


d. 


(30) 


where 

f ( • )  =  i  I  r  r  +  •  (z  +  d)  +  : l/kR 
s<  ■  >  - 

R  =  Vr^  +  ( z  +  d )  2 
=  kR. 
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The  parameter  p  is  introduced  so  that  we  have  an  adjustment  on  the  rate  of 
descent  of  the  function  away  from  the  saddle  point.  When  p  is  large,  the 
magnitude  will  decrease  quickly  on  either  side  of  the  saddle  and  our  approxi¬ 
mation  will  be  good.  The  stationary  points  of  the  exponent  are  given  by, 


or 


or 


i 

kR 


[r  + 


3d. 

3> 


^  (z  +  d)]  =  0 


r 


(z  +  d) 


(31) 


At  this  point,  let  us  examine  the  geometry  of  the  reflection  process. 

In  figure  5,  we  have  shown  the  ray  path  {the  normal  to  the  wavefront)  for  a 
wave  with  horizontal  wave  number  y0.  The  length  of  the  projection  of  the 
legs  of  the  ray  path  onto  the  interface  is  y<j(z  +  d )/ e -  Thus,  t.ie  displacement 
of  the  field  on  reflection  is  given  by  the  difference. 


r  -  ^  (z  +  d)  =  6 


which  is,  by  equation  (31), 


=  _  jJL 
37 


(32) 


This  quantity  is  usually  known  as  beam  displacement  because  of  its  derivation 
in  terms  of  bounded  beams^,  but  any  nonplanar  wavefront  exhibits  a  similar 
displacement  so  field  displacement  may  be  a  more  appropriate  term. 

Notice  that,  if  we  had  not  included  the  phase  of  the  reflection  coefficient 
in  the  exponent,  we  would  have  obtained  the  naive  result  that  there  is  no 
displacement  of  the  field  on  reflection.  The  most  obvious  case  in  which  there 
is  a  displacement  of  the  field  is  a  reflector  that  is  a  positive  gradient 
half-space.  A  ray  that  enters  the  gradient  half-space  will  eventually  exit 
further  downrange  as  it  is  refracted  upwards  by  the  sound  speed  gradient. 

The  displacement  given  by  equation  (32)  is  identical  to  the  displacement  given 
by  ray  theory.  There  is  also  a  displacement  on  reflection  from  a  homogeneous 
half-space.  This  is  a  wave  phenomeon  as  ray  theory  predicts  no  such  displace¬ 
ment.  This  displacement  has,  however,  been  observed  experimentally. 3 
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By  the  method  of  steepest  descent  the  integral's  value  resulting  from 
one  saddle  point  is, 


,°“eCf<'o)/3i^G<,°)  (33) 

where  only  the  first  term  of  the  series  expansion  of  G  about  the  saddle  point 
has  been  retained.  As  long  as  kR»l,  this  will  be  a  valid  approximation. 
Substituting  equation  (33)  for  the  integral  in  equation  (27),  we  can  write 
the  field  as. 
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"V 


+  e(z  +  d)  -wt] 
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+  d)  +  6] 


(34) 


The  significance  of  the  factor  under  the  radical  is  not  immediately  obvious  but 
it  represents  the  geometrical  spreading  loss  over  the  propagation  path  including 
the  displacement  on  reflection.  This  can  be  seen  in  the  construction  shown 
in  figure  6. 

The  geometrical  spreading  loss  is  given  by  the  square  root  of  the  ratio  of  the 
elemental  areas  m  and  e0m0: 


rB2  17  [s  (2  +  d)  +  6] 

(The  field  intensity  is  inversely  proportional  to  the  cross-sectional  area  of 
the  ray  bundle,  but  ;  is  field  amplitude  which  is  proportional  to  the  square 
root  of  field  intensity.)  Consequently,  the  first  level  of  approximation 
(equation  (34))  by  the  method  of  steepest  descent  represents  the  field  as 
ordinary  plane  wave  reflection  with  an  appropriate  displacement  on  reflection 
and  geometrical  spreading.  If  the  reflection  coefficient  is  considered  to  be 
a  constant,  then  there  will  be  no  field  displacement  and  equation  (34)  reduces 
to  the  simple  plane  wave  result, 

i  [iQr  +  :  ( z  +  d)  -  u>t] 


Vr2  +  (z  +  d)2  ' 

So  far  in  this  section,  we  have  established  the  relationship  between  the  full 
wave  solution  for  reflection  and  the  ray  theory  solution  (geometrical  spreading) 
with  plane  wave  reflection.  The  general  criterion  for  validity  of  the  latter 
technique  is  that  the  quantity  kR  be  much  greater  than  unity.  Another  way  of 
considering  this  criterion  is  to  define  a  local  wavelength  angle  .  At  the 
range  in  question  (the  range  to  the  reflector,  for  example),  this  is  the  angle 
subtended  by  an  arc  one  wavelength  long  along  the  spherical  wavefront.  Thus, 
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kR  »  1 


9x  «  2n 


becomes , 


where 


Therefore,  the  validity  of  the  plane  wave  assumption  depends  not  on  the  wave- 
front  curvature  alone  but  on  the  curvature-wavelength  product.  If  one  wavelength 
of  arc  along  the  incident  spherical  wave  represents  at  least  a  substantial  part 
of  a  whole  circle,  then  spherical  wave  theory  should  be  considered. 

As  we  have  noted,  this  is  an  area  in  which  some  caution  must  be  exercised. 
While  we  have  examined  the  full  solution  and  an  intermediate  approximation 
in  which  the  more  pronounced  effects  of  a  nonplanar  incident  wavefront  are 
considered,  we  must  not  forget  the  underlying  assumptions.  For  example,  we 
are  depending  on  slow  variation  of  the  magnitude  of  the  reflection  coefficient. 

We  are  also  neglecting  the  magnitude  variation  of  the  /y/e  factor  in  the 
integrand.  Both  of  these  assumptions  are  violated  in  certain  regions  (such 
as  near  a  critical  angle)  and  we  must  allow  for  these  violations.  In  the 
next  section,  we  will  consider  an  example  in  which  we  can  see  the  effects  of 
some  of  these  problem  regions  and  also  how  to  extract  meaningful  solutions 
in  spite  of  the  difficulties. 


REFLECTION  FROM  A  THIN  LAYER 


To  illustrate  the  effects  of  poles  and  a  branch  line  integral  while 
maintaining  some  simplicity  in  mathematics,  let  us  consider  reflection  from  a 
finite  fluid  layer  over  a  rigid  reflector  as  pictured  in  figure  7.  The  plane 
wave  reflection  coefficient  for  this  case  (evaluated  at  z=0)  is, 

3  cos  3,h  +  ibe,  sin  ts,h 

_  O  1  i _ 1  /• 3C 


3  cos  B,h  -  ibe,  sin  t-.h 
o  1  11 


where 


=  'o/cl 


*o,l  *  V(“Z/co.1>  -  '-2  • 

For  real  angles  of  incidence  (yWc0),  the  reflection  coefficient  magnitude  is 
always  one  since  all  the  incident  energy  is  returned  either  by  total  reflection 
at  z=0  or  by  reflection  from  the  rigid  surface  at  z=-h.  The  reflected  field  is 
given  by  equations  (27  and  (35).  Since  equation  (35)  is  symmetric  in  but 
not  in  30,  there  will  only  be  one  branch  line  integral  and  that  will  be 
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associated  with  the  pole  at  6o=0.  In  addition,  we  can  expect  a  number  of  poles 
for  complex  values  of  y  that  are  associated  with  modes  partially  trapped  in 
the  finite  layer.  Since  ci>c0,  there  will  be  no  completely  trapped  modes. 

Before  we  examine  the  wave  theory  solution,  let  us  consider  what  ray 
theory  would  predict.  There  should  be  one  path  directly  from  the  source  to 
the  upper  surface  of  the  layer  and  back  to  the  receiver.  Also,  there  will  be 
an  infinite  number  of  paths  that  reflect  from  both  the  rigid  reflector  and 
from  the  underside  of  the  layer's  upper  surface  some  number  of  times  before 
exiting  the  layer.  These  paths  will  be  progressi vely  less  important  as  the 
number  of  internal  bounces  increases  since  some  energy  is  lost  at  each  reflection 
with  the  upper  surface  of  the  layer.  While  we  should  expect  some  differences 
in  the  wave  theory  treatment,  we  would  not  expect  the  descent  method  to  predict 
only  a  single  path,  reflected  and  displaced. 

As  a  starting  point  for  any  of  the  wave  solutions,  figure  8  shows  two 
possible  configurations  of  the  integrand  of  equation  (27)  in  the  complex  wave 
number  (y)  plane.  In  the  first  case,  deformation  of  the  original  integration 
contour  leads  to  a  branch  line  integral  only,  while  in  the  second  case, 
residues  from  each  of  tiie  poles  must  be  added  to  a  somewhat  simpler  branch 
line  integral.  These  techniques  are  straightforward;  however,  we  must  specify 
the  particular  branch  of  s0.  In  each  of  the  regions  (of  the  right  half-plane) 
of  figure  8,  the  branches  of  6q  are  as  follows, 


\ 

) 

6o  =  ■  1 

II) 

=  +1 

III) 

6o  =  +  1 

k--2' 

For  convenience  in  the  following  discussion,  we  will  refer  to  the  sheet  of  the 
integrand  on  which  the  positive  root  of  3o  is  used  as  the  upper  sheet  and  the 
sheet  on  which  the  negative  root  is  used  as  the  lower  sheet.  By  this  definition, 
the  poles  shown  in  figure  8  are  on  the  upper  sheet.  Also,  this  definition  implies 
a  branch  cut  extending  along  the  real  y-axis  to  the  right  of  the  pole  at  eo=0- 
This  is  the  conventional  branch  used  by  computer  complex  square  root  routines. 

We  will  now  consider  a  solution  by  means  of  the  method  of  steepest  descent. 
Although  the  numbers  are  themselves  unimportant,  we  will  use  the  following  case 
to  illustrate  the  method: 

cQ  =  1500  m/s 

Cj  =  2000  m/s 

h  =  100  m 
z  -  d  =  50  m 
cj/o2  =  0.5 
f  =50  Hz. 
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The  first  step  involves  locating  the  saddle  points  or  stationary  phase  points 
by  means  of  equation  (31).  Since  the  right-hand-side  of  this  equation  is 
independent  of  the  source-receiver  geometry,  we  can  plot  this  function 
separately  and  then  plot  the  left-hand-side  parametrically  in  range.  The  points 
of  intersection  determine  the  stationary  phase  points.  Such  a  plot  is  shown 
in  figure  9.  Most  of  the  time,  each  pole  (evidenced  by  the  peaks  in  the  field 
displacement  curve)  gives  rise  to  a  pair  of  stationary  phase  points  and  then 
there  is  one  additional  stationary  point  to  the  right  of  all  of  the  pole  pairs. 
Thus,  the  influence  of  each  pole  (the  partially-trapped  modes)  is  felt  through 
a  pair  of  saddle  points  and  the  path  related  to  the  reflection  from  the  top  of 
the  layer  is  accounted  for  by  the  remaining  stationary  point. 

Once  the  stationary  phase  points  have  been  located,  the  contours  of  steepest 
descent  must  be  found.  In  case  it  is  necessary  to  cross  poles  or  branch  cuts 
in  order  to  deform  the  original  contour  into  the  descent  path,  the  appropriate 
residues  or  BLIs  must  be  added  to  the  saddle  point  solutions.  The  descent 
paths  for  our  example  (and  r=300  m)  are  shown  in  figure  10  for  only  the  reflec¬ 
tion  coefficient  and  exponential  factors  of  equation  (27).  The  branch  cut  of 
figure  8b  is  used  because  the  descent  contour  must  be  continuous.  We  can  see 
that  in  the  process  of  deforming  the  original  contour  (the  real  y-axis)  into 
the  descent  path,  no  poles  are  crossed  and  the  branch  cut  is  avoided.  Hence, 
the  solution  is  complete  with  only  the  saddle  points  considered.  Also,  notice 
the  path  of  descent:  between  each  pair  of  saddle  points  associated  with  a 
pole,  the  path  dips  into  a  zero  rather  than  going  to  negative  infinity.  As 
can  be  seen  from  figure  9,  as  the  range  increases  the  two  stationary  phase 
points  of  a  pair  draw  closer  together.  Consequently,  the  descent  from  either 
saddle  into  the  zero  becomes  less  steep  and  more  terms  may  be  necessary  in  the 
approximate  solution.  Fortunately,  the  modes  which  cause  this  trouble  first 
are  the  modes  that  correspond  to  steeper  angles  of  penetration  into  and  out  of 
the  layer  and,  therefore,  lose  more  energy  than  the  "shallower"  modes. 

As  we  have  mentioned,  the  assumption  that  the  /y/s  factor  is  slowly  varying 
is  not  always  a  good  assumption.  Figure  11  shows  the  descent  paths  and  saddle 
points  for  the  entire  integrand  of  equation  (27).  While  the  pattern  is  very 
similar  to  that  in  figure  10,  there  are  some  important  differences.  First, 
the  saddle  points  have  moved  slightly  (both  figures  10  and  11  are  exaggerated 
in  the  vertical  direction  by  a  factor  of  ten).  Since  the  /y/e  factor  could 
be  expressed  as  a  complex  exponential,  the  resulting  exponent  will  add  another 
term  to  the  stationary  phase  point  relation,  equation  (31),  which  would  predict 
the  shift.  Second,  as  a  consequence  of  the  saddle  point  shift  associated  with 
the  most  lossy  mode  (farthest  to  the  left  and  away  from  the  real  axis),  the 
descent  path  has  made  a  definite  change.  The  path  passes  through  only  one 
of  the  saddle  points  now  and  does  not  venture  into  the  zero  below  this  pole. 

From  the  stationary  phase  plot  in  figure  9,  we  would  expect  that  as  the 
range  increases,  first,  there  would  be  two  saddle  points  for  a  given  mode, 
then  at  some  particular  range,  they  would  somehow  coalesce  into  a  single  saddle 
which  would  then  disappear  for  greater  ranges.  This  is  not,  in  fact,  the  case 
as  we  see  in  the  sequence  in  figure  12.  Here,  we  have  plotted  the  contours  of 
the  entire  integrand  for  three  ranges  -  200,  300,  and  400  m  -  in  the  vicinity 
of  the  left-most  pole-zero  combination.  At  200  m,  the  two  saddles  are  nearly 
symmetrical  and  the  descent  path  passes  through  both  saddles  and  the  zero. 

At  300  m,  the  >  ./  factor  has  warped  the  contours  so  that,  because  the  saddles 
have  drawn  closer  together,  the  steepest  descent  path  swings  positively  upward 
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in  either  direction  away  from  the  right-hand  saddle.  Thus,  the  left-hand 
saddle  is  bypassed  in  the  contour  deformation.  As  the  range  increases  further, 
the  right-hand  saddle  turns  so  that  the  path  becomes  horizontal  while  the  left- 
hand  saddle  turns  toward  the  vertical.  The  contour  plot  at  r  =  400m  illustrates 
this  state  and  shows  that,  even  though  no  stationary  phase  point  is  predicted 
by  equation  (31)  through  figure  8,  one  does,  in  fact,  exist  with  a  related 
descent  path.  This  illustrates  a  short-coming  in  equation  (31).  The  steepest 
descent  method  will  produce  the  correct  answers  only  if  the  relevant  saddle 
points  are  properly  identified.  Equation  (31)  gives  good  approximations  for 
most  of  these  saddle  points  but,  in  some  instances,  one  or  two  may  be  repre¬ 
sented  improperly  or  missed  entirely. 

Pernaps  it  is  safer  practice  to  solve  for  the  field  by  a  hybrid  techniaue 
in  which  we  deform  the  contour  past  all  of  the  poles  but  still  let  it  pass 
through  the  right-most  saddle  point  in  order  to  avoid  the  branch  line.  In 
this  way,  the  solution  becomes  a  sum  of  residues  and  a  series  approximation  at 
one  saddle  point. 

This  solution  process  tells  us  something  about  the  resultant  signal's 
phase.  Since  this  modeling  effort  was  undertaken  to  develop  a  tool  useful 
for  nonlinear  signal  processing,!  we  are  particularly  interested  in  the  phase 
of  the  field.  From  figure  9,  we  can  see  that  the  phase  change  is  rapid  in  the 
vicinity  of  the  poles  of  the  reflection  coefficient.  This  will  lead  to 
problems  in  phase  reconstruction  of  sampled  signals  unless  the  direction  of  the 
phase  change  is  known.  In  order  to  have  a  positive  field  displacement  (negative 
displacements  have  not  been  observed  in  any  of  the  test  cases),  the  rate  of 
change  of  phase  with  respect  to  wave  number  must  be  negative  by  equation  (32). 
Since  horizontal  wave  number  is  directly  proportional  to  frequency  at  constant 
ray  geometry,  the  rate  of  change  of  phase  with  respect  to  frequency  must  also 
be  negative: 


where 


cv  =  ray  equivalent  vertex  velocity. 

If  the  field  displacement  is  always  positive  in  real  physical  reflection 
problems,  then  we  know  that  the  phase  progression  must  be  characterized  by 
rapid  decreases  in  the  vicinity  of  the  reflection  coefficient  poles.  This 
principle  would  be  of  great  value  in  phase  reconstruction  of  noisy,  sampled 
signals.  See  reference  1  for  a  more  complete  discussion  of  phase  reconstruction 
problems  and  techniques. 

To  summarize  this  section,  let  us  review  the  levels  of  approximation  we 
now  have  avai lable: 

1.  By  neglecting  the  rate  of  phase  change  of  the  reflection  coefficient, 
we  locate  a  single  stationary  phase  point  which  corresponds  to  specular  reflec¬ 
tion  with  no  displacement.  Any  multipath  phenomena  associated  with  the  reflector 
are  ignored. 
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2.  If  we  locate  all  of  the  stationary  phase  points  (saddle  points)  that 
result  from  considering  the  reflection  coefficient  phase  variation,  then  we 
can  write  approximations  for  the  multiple  paths  of  propagation  from  the 
reflector.  The  possibility  does  exist  for  misinterpreting  one  or  more  of 
these  points  because  of  the  effects  of  other  factors  in  the  integrand. 

3.  We  can  also  consider  the  entire  integrand  in  order  to  fix  the  location 
of  the  saddle  points  and  descent  path.  This  will  improve  the  accuracy  of  the 
saddle  point  approximations  at  the  expense  of  substantially  increased  computa¬ 
tional  labor. 

4.  In  some  cases,  the  poles  themselves  may  be  relatively  straightforward 
to  locate,  in  which  case  many  of  the  saddle  point  approximations  can  be  replaced 
by  residues. 

5.  Finally,  the  integral  can  be  numerically  integrated  as  it  appears  in 
equation  (27).  This  procedure  can  be  expensive  particularly  when  the  impulse 
response  of  the  reflector  is  desired  since  an  additional  integration  over 
frequency  is  required. 


NUMERICAL  IMPLEMENTATION 


Because  most  of  the  procedures  discussed  above  are,  in  practice,  applied 
numerically  rather  than  analytically,  we  will  now  consider  some  numerical 
techniques  for  solving  the  field  integral.  In  order  to  start  with  as  close 
to  an  exact  solution  as  possible,  an  algorithm  has  been  developed  based  on 
Bucker's  technique8  to  integrate  the  field  integral,  equation  (27)  directly. 

Since  variations  in  the  integrand  are  irregular  in  that  along  some  portions 
of  the  contour  the  function  is  smooth  and  along  other  portions  the  function 
changes  rapidly,  the  integration  procedure  should  be  adaptive.  This  is  accom¬ 
plished  by  integrating  over  a  segment  (in  this  case,  by  Gauss-Legendre 
Quadrature  with  n=5)  and  comparing  this  result  with  the  sum  of  two  integrations 
over  the  two  halves  of  the  same  segment.  If  the  difference  is  negligible,  the 
step  size  is  increased  for  the  next  segment.  Otherwise,  the  step  size  is 
halved  and  the  integrations  are  performed  again.  In  this  way,  the  step  size 
is  continually  adjusted  according  to  the  rate  of  variation  of  the  integrand. 


As  in  some  of  the  processes  to  follow,  a  small  amount  of  attenuation  is 
introduced  so  as  to  effectively  move  the  integration  contour  slightly  off  the 
real  axis.  This  avoids  the  complications  associated  with  integration  through 
the  pole  at  ;-o=0.  The  attenuation  is  treated  as  viscous  attenuation  (constant 
loss  per  cycle)  by  adding  a  small  imaginary  part  to  the  wave  number.  Actually, 
this  is  done  by  making  the  sound  speed  complex  in  such  a  way  as  to  prescribe 
the  appropriate  wave  number. 


Besides  direct  integration  of  equation  (27),  another  solution  technique 
involves  applying  the  FFT  algorithm  to  the  Fourier  transform  integral  of 
equation  (28).  This  approach  is  particularly  convenient  if  we  would  like  the 
field  at  many  different  ranges. 9  Direct  numerical  integration  only  gives  one 
range  value  per  integration  while,  for  example,  a  512  point  FFT  would  produce 
field  values  at  512  evenly  spaced  ranges.  The  principal  disadvantage  of  this 
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method  is  that  the  integrand  (or,  rather,  the  F(y)  function  of  equation  (28)) 
varies  quite  rapidly  over  some  values  of  y  thereby  requiring,  as  a  consequence 
of  the  sampling  theorems,  fine  sampling  of  the  entire  function.  Much  of  this 
computational  work  is  wasted  because  the  function  can  be  smooth  over  other, 
large  segments  of  y  where  the  sampling  does  not  need  to  be  as  frequent.  Since 
the  sampling  must  be  uniform  over  the  entire  range  of  r  (a  restriction  imposed 
by  the  FFT),  the  regions  of  rapid  variation  determine  the  overall  sampling  rate 
and,  therefore,  the  transform  size.  Without  any  special  treatment  then,  this 
method  can  involve  expensive  computation  of  many  more  points  than  are  necessary 
to  define  the  function  and  very  large  transform  sizes. 

It  is  possible,  however,  to  use  digital  filter  theory  to  greatly  increase 
the  efficienty  of  the  FFT  technique.  Let  us  consider  the  development  of  such 
a  process.  Equation  (28)  can  be  written  in  a  slightly  different  form  so  that 
the  integrand  function  F  and  the  field  function  are  a  Fourier  transform  pair. 


?(r)  =  /  F(Y)eiYrd-,  (36) 


Here,  the  function  F  can  be  determined  from  equation  (27). 

Since  and  F  are  a  Fourier  transform  pair,  we  can  filter  either  function 
in  either  the  wave  number  or  the  range  domain  and  also  we  can  expect  that  once 
the  functions  are  sampled,  the  theorem,  relating  sampling  rate  and  bandwidth, 
will  apply.  We  shall  see  later  that  this  is  true.  Meanwhile,  it  is  sufficient 
to  say  that  if  we  can  restrict  the  "bandwidth"  of  :  in  the  range  domain  (that  is, 
restrict  the  span  of  relevant  ranges),  then  we  can  sample  the  function  F  at  a 
much  lower  rate  than  would  normally  be  required  to  properly  define  the  function. 
If  we  knew  p,  we  could  restrict  its  bandwidth  by  setting  all  values  of  ;  above 
some  maximum  range  R  to  zero.  This  is  equivalent  to  filtering  F(y)  in  such  a 
way  that  the  high  rate  variations  are  suppressed  (low  pass  filtering).  To 
perform  this  operation,  then,  we  must  locally  sample  F(y)  densely  enough  to 
describe  the  function  variations,  filter  these  sections,  and  reduce  the  sampling 
rate  according  to  the  sampling  relations  derived  below.  We  then  would  have  an 
efficient  technique  for  field  integral  solution  by  FFT. 

Now  let  us  examine  the  conversion  of  equation  (36)  into  a  discrete  Fourier 
transform  suitable  for  the  FFT  algorithm.  We  know  that  the  FFT  operates  on  a 
finite  length  series  of  samples  of  a  function.  Since  a  periodic  function  can 
be  exactly  represented  by  a  Fourier  series  expansion  in  place  of  a  continuous 
Fourier  transform,  we  can  convert  the  integral  of  equation  (36)  into  a  summation 
if  we  are  willing  to  accept  a  field  function  periodic  in  range.  If  we  select  a 
period  R  that  is  greater  than  the  maximum  range  of  interest  this  will  not 
present  a  problem  in  itself.  We  will  have  to  filter  F(.)  as  discussed  above 
so  that  :(r)  is  zero  for  r'-R  prior  to  sampling  F(y).  If  we  do  not  filter, 
then  (r)  will  be  distorted  by  aliasing. 
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To  make  the  summation  finite,  we  can  taper  F(y)  to  zero  beyond  some 
maximum  r.  This  will  cause  smoothing  of  the  function  $  (i.e.,  low  pass 
filtering)  in  range  but  this  is  actually  desireable  in  a  wave  theory  solution. 
At  this  point,  we  have  reduced  equation  (36)  to  the  following  form, 


N-l 

<t>  (k&r)  = 

n=0 


F(nAY)e-iknAr^  ay 


(37) 


where 

k,n  =  integer  indexes 
Ar  =  sample  interval  of  range 
Ay  =  sample  interval  of  wave  number 
N  =  transform  size  such  that  NAr  =  R. 

One  more  modification  is  now  necessary  before  the  FFT  can  be  applied.  Let, 

kr-ArAy  =  2^kn/N 


or 


ArAy  =  2it/N 


so  that. 


N-l 

4>(kAr)  =  Ay  ^  ^  F(nAy)e  i2nkn/,N 
n=0 


(38) 


=  Ay  FFT  (F(nAy) ]  (39) 

Equation  (38)  is  essentially  a  statement  of  the  sampling  restrictions  for  this 
problem.  Since  we  can  determine  the  maximum  range  R  and 

R  =  NAr 

then  the  increment  Ay  is  no  longer  arbitrary  but,  by  equation  (38), 

Ay  =  2-t/R  (40) 

and  therefore,  the  relevant  span  of  wave  number  is, 

r  =  Nay  =  2  r/Ar  (41) 
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In  practice,  we  would  select  the  transform  size  and  either  R  or  Ar.  This 
fixes  the  rest  of  the  parameters  and  also  means  that  only  as  many  F(y)  points 
are  needed  as  range  points  are  desired.  Once  again,  the  function  F(y)  must 
be  initially  sampled  in  such  a  way  that  its  variations  are  locally  well  defined 
close  sampling  in  regions  of  rapid  variation  and  infrequent  sampling  in  regions 
of  slow  variation  -  but  then  the  function  is  filtered  and  resampled  uniformly 
and  at  a  low  rate  given  by  equation  (40).  The  simplest  filter  to  use  is  a 
convolutional  filter  designed  by  the  window  method  based  on  a  function  in  the 
range  domain  that  is  unity  between  r=Q  and  r=R  and  zero  otherwise.  The 
transform  of  this  range  function  is  truncated  by  one  of  the  standard  windows 
in  the  wave  number  domain  and  then  applied  in  a  point-by-point  convolution  to 
the  F(y)  function.  The  convolution  is  only  performed  at  those  locations  that 
will  be  new  sample  points  of  F(y). 

We  have  seen  how  the  FFT  can  be  used  to  compute  many  range  points  of  the 
field  function  efficiently.  Since  this  study  is  primarily  concerned  with  time- 
domain  results,  we  will  now  examine  the  use  of  the  FFT  in  obtaining  time  series 
waveforms  from  the  field  integral  solutions.  Before  we  considered  any  of  the 
field  integral  solutions,  we  reduced  the  wave  equation  to  a  Helmholtz  equation 
by  assuming  a  time  dependence  of  e-iwt.  This  reduction  could  also  have  been 
done  by  taking  the  Fourier  transform  of  equation  (1).  We  can  thus  see  that, 
in  order  to  obtain  the  time  series,  we  need  only  take  the  inverse  Fourier 
transform  of  the  field  function  Hence, 


co 


-CO 


where 


p  =  pressure  as  a  function  of  time. 

This  is  basically  the  same  result  as  equation  (18).  Notice  that  the  independent 
variable  in  the  function  t>  is  frequency  w  rather  than  range.  As  a  result,  we 
must  compute  the  field  function  at  many  frequencies  and  a  single  range  for  each 
p(t)  function. 


As  before,  we  can  recast  equation  (42)  in  a  form  suitable  for  using  the  FFT 


p(nAt)  =  Af 


N-l 

L 

k=0 


-  i  2  -rkn/N 


=  :  f  F FT- 1  [ :  (k.\f)  ] 


(43) 
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where 


Af  =  discrete  form  of  dw/2u 
N  =  l/(AtAf). 

The  pressure  response  p  is  the  impulse  response  whose  spectrum  is  given  by  4>(u>). 
The  implied  impulse  strength  can  be  found  by  evaluating  equation  (43)  at  n  =  0 
(which  is  to  say  t=0), 

P0  ■  NAf 


or 


PQAt  =  NAfAt  =  1 

which  is  the  impulse  strength.  We  will  find  it  convenient  to  make  all  of  the 
time-domain  (impulse  response)  calculations  in  terms  of  impulse  strength  so 
that  the  results  are  independent  of  transform  size.  We  will  then  compute, 

p(nAt)At  =  -jj-FFT'^kAf)]  (44) 

Before  we  leave  our  discussion  of  numerical  techniques,  let  us  consider 
two  more  aspects  of  the  FFT  techniques.  One  of  these  we  have  mentioned  before 
and  that  is  the  inclusion  of  the  source  waveform  by  spectrum  multiplication. 

In  particular,  we  will  be  modeling  a  source  in  the  next  section  that  can  be 
approximated  by  an  exponential  function, 

s ( t )  =  Ae_at  (t  >  0) 


where 


a  =  decay  time  constant. 

First,  we  would  like  to  normalize  the  source  so  that,  when  applied  by  convolu¬ 
tion  to  a  uniform  level,  it  does  not  change  the  power  of  the  time  function. 
This  can  be  assured  by  letting, 


1  = 


/ 


Ae'atdt 


or 

A  =  x  . 

Therefore,  we  will  use  a  normalized  source  function  of. 


s ( t )  =  <e"lt  (t  >  0) 


(45) 
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although  it  will  be  applied  by  multiplying  the  field  function  $(w),  by  the 
source  spectrum  S(w), 


S(«u) 


-at  ioot 
ae  e 


dt 


and  the  desired  pressure  response  is  then, 

■1/ 


r(t)At  =  jJ-FFT' 


’  +  4)  J 


where 


(46) 


r  =  simulated  received  signal. 

Finally,  the  truncation  of  t>(w)  beyond  some  maximum  frequency  introduces 
spurious  oscillations  into  the  p(t)  function.  The  truncation  is  necessary 
because  of  the  finite  length  of  the  transform;  however,  the  truncation  does 
not  have  to  be  abrupt.  If  we  taper  the  upper  edge  of  the  $  spectrum  by  some 
function  (perhaps  half  a  cycle  of  a  raised  cosine)  prior  to  transforming,  we 
can  reduce  these  oscillations  substantially.  In  practice  then,  the  simulated 
received  signal  is  computed  as  follows. 


>’(  t). 


t  =  — 

1  N 


WU)*(u) 


(47) 


where 


W(j)  =  weighting  function. 


APPLICATIONS 


One  of  the  first  decis.ons  that  must  be  made  in  using  these  models  for 
transient  propagation  is  whether  to  use  a  model  based  on  plane  waves  or 
spherical  waves.  Predictions  based  on  plane  wave  theory  are  generally  fast 
and  easily  applied  to  complicated  multipath  propagation.  On  the  other  hand, 
there  are  instances  in  which  the  plane  wave  assumptions  lead  to  grossly 
incomplete  results.  In  this  section,  we  will  duscuss  those  effects  predicted 
by  both  the  plane  and  spherical  wave  theories,  develop  the  criteria  for  decision 
between  the  two,  and  illustrate  prediction  by  the  full  wave  theory. 
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Generally,  the  plane  wave  theory  is  used  as  an  adjunct  to  ray  or  geometrical 
acoustics.  Spreading  losses  and  refraction  in  the  water  are  accounted  for  by 
ray  tracing  and  the  ray  path  is  modified  in  magnitude  and  phase  at  the  ocean 
floor  according  to  the  plane  wave  reflection  coefficient.  We  will  add  to  this 
the  possibility  of  a  displacement  (or  even  multiple  displacements)  on  reflection 
from  the  ocean  floor  because  this  phenomenon  is  included  by  the  degenerate  form 
of  the  full  solution,  equation  (34).  This  displacement  is  a  legitimate  feature 
of  the  field  as  long  as  the  magnitude  of  the  reflection  coefficient  is  slowly 
varying  with  wave  number. 

Plane  wave  theory  then  includes  the  following  effects:  time  delay, 
refraction,  and  spreading  in  the  water  by  ray  acoustics;  phase  and  magnitude 
changes  and  field  displacement  on  reflection  from  the  ocean  floor;  and  phase 
distortion  at  ray  turning  points.  The  last  effect  amounts  to  a  uniform  -r/2 
phase  shift  at  a  turning  point  (caused  by  refraction)  and  is  evidenced  by 
formation  of  doublets  in  the  time  series  (a  single  positive  impulse  becomes  a 
positive  impulse  followed  immediately  by  a  negative  impulse).  This  doublet 
formation  is  obvious  in  the  time  series  and  can  be  used  to  identify  refraction 
turning  points  either  in  the  water  column  or  in  the  ocean  floor.  Host  of  the 
other  effects  are  more  difficult  to  use  in  direct  interpretation  of  measure¬ 
ment  signals  but  the  time  delay  is  valuable  for  isolating  distinctly  different 
propagation  paths. 

Figure  4  illustrates  the  construction  of  a  plane  wave  model  that  treats 
a  number  of  different  paths.  Each  path  is  characterized  by  a  water  path  delay 
and  loss,  and  a  bottom  reflection  characterized  by  the  plane  wave  reflection 
coefficient.  Since  the  paths  are  not  mutually  dependent,  they  can  be  added  or 
deleted  as  necessary  to  represent  any  aspect  of  an  ocean  environment.  This 
feature,  coupled  with  the  dependence  of  the  reflection  coefficient  only  on 
angle,  yields  a  model  that  can  be  easily  adapted  to  different  situations. 

This  approach  should  be  used  whenever  it  is  applicable  because  it  is  not  only 
flexible  but  also  very  inexpensive  compared  to  the  spherical  wave  solution. 

As  we  will  outline  more  specifically,  there  are  important  applications  in 
which  the  plane  wave  solution  is  inadequate.  Besides  the  effects  covered  by  the 
plane  wave  model,  the  spherical  wave  solution  includes  boundary  waves  of  several 
types  and  ducted  propagation  in  layered  or  refracting  bottoms.  Particularly 
at  low  frequency  or  in  shallow  water,  these  mechanisms  can  contribute  a  sub¬ 
stantial,  if  not  dominant,  portions  of  the  received  energy. 

In  the  case  in  which  the  ocean  floor  is  modeled  as  a  homogeneous  fluid, 
an  interface  wave  can  be  excited  which  travels  horizontally  while  continuously 
radiating  energy  back  into  the  overlying  ocean.  This  interface  wave  is  excited 
by  energy  incident  to  the  boundary  at  the  critical  angle  and  the  re-radiation 
appears  to  be  coming  from  the  boundary  at  that  same  critical  angle.  Because 
of  the  continuous  re-radiation  associated  with  this  boundary  wave,  the  energy 
in  the  advancing  wave  decays  much  more  rapidly  than  the  specularly  reflected 
wave.  (The  amplitude  of  the  interface  wave  drops  at  approximately  r“2  while 
that  of  the  reflected  wave  drops  at  r- 1.)  As  a  result,  this  interface  wave 
is  only  significant  at  shorter  ranges,  but  it  can  change  the  propagation 
characteri sties  substantially  at  ranges  where  the  specular  reflection  is  close 
to  the  critical  angle  (figure  13).  This  type  of  wave  is  also  known  a 
critically  refracted  wave  or  lateral  wave. 
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If  the  fluid  representing  the  ocean  floor  is  layered  or  has  a  sound 
speed  profile  that  produces  a  duct,  then  partially  trapped  or  ducted  propaga¬ 
tion  can  take  place.  Ray  theory  can  be  used  to  model  ducting  but,  especially 
when  the  acoustic  wavelength  is  of  the  order  of  the  duct  width,  the  integral 
wave  theory  provides  much  better  representation  of  the  ducted  energy.  Typically, 
there  will  be  a  number  of  stationary  phase  points  corresponding  to  the  poles 
which,  in  turn,  represent  the  leaky  modes  in  the  duct. 

Additional  types  of  waves  appear  if  the  bottom  model  is  made  of  an  elastic 
material  (that  is,  a  solid).  At  least  two  other  types  of  interface  waves  can 
propagate  in  this  case:  a  shear  wave  analogous  to  the  compressional  interface 
wave  for  the  fluid  reflector  and  a  Stoneley  wave  in  which  the  particle  motion 
is  confined  to  the  vertical  plane  containing  the  source  and  receiver.  The 
Stoneley  wave  has  its  own  pole  in  the  reflection  coefficient  for  some  complex 
wave  number  and  it  loses  energy  (as  does  the  other  shear  interface  wave)  at  a 
rate  similar  to  that  of  the  compressional  interface  wave. 

Ducts  caused  by  solid  layers  or  gradients  in  the  shear  sound  speed  give 
rise  to  ducted  shear  propagation  or  Love  waves.  One  noteable  characteristic 
of  either  shear  or  compressional  ducted  propagation  is  that  it  is  dispersive: 
different  frequencies  travel  with  different  speeds-  This  results  in  a  dis¬ 
tortion  of  the  time  waveshape  that  the  plane  wave  technique  is  unable  to 
describe.  This  also  can  make  replica  deconvolution  (inverse  filtering) 
difficult  since  the  dispersion  effectively  distorts  the  source  waveform. 

An  important  aspect  of  shear  propagation  in  real  ocean  sediments  is 
attenuation.  Shear  attenuation  is  generally  much  greater  than  compressional 
attenuation  because  sediments  are  much  closer  to  fluids  than  solids  even  though 
they  do  have  some  rigidity.  Consequently,  if  solid  properties  are  accounted 
for  by  the  model,  attenuation  must  be  included.  Otherwise,  the  shear  propa¬ 
gation  can  be  unrealistically  good. 

Having  summarized  the  scope  of  both  the  plane  wave  and  the  spherical 
wave  solutions,  let  us  now  outline  the  criteria  for  deciding  which  to  use. 

As  we  have  shown,  the  mathematical  criteria  for  validity  of  the  plane  wave 
approximation  are  yr>>l  (so  that  the  asymptotic  form  of  the  Hankel  function 
can  be  used)  and  kR>>l.  These  are  roughly  equivalent  and  the  latter  criterion 
amounts  to  insuring  that  a  wavelength  is  small  compared  to  the  wavefront 
circumference  at  the  radius  of  interest.  Another  way  of  saying  this  is  that 
2-t  times  the  wavefront's  radius  of  curvature  must  be  much  greater  than  the 
acoustic  wavelength. 

If  we  desire  to  simplify  the  plane  wave  theory  even  further  by  neglecting 
field  displacement  and  multiple  stationary  points,  we  need  some  additional 
restrictions.  First,  trapped  propagation  in  the  bottom  must  be  negligible. 
Trapped  modes  are  evidenced  by  poles  in  the  reflection  coefficient  and  these 
poles  will  either  lead  to  multiple  stationary  phase  points  or  residues. 

Second,  the  geometry  must  be  such  that  the  critically  refracted  or  Stonely 
interface  waves  are  not  favorably  excited.  In  order  to  avoid  exciting  the 
critically  refracted  wave,  the  geometry  must  be  such  that  the  ordinary 
reflected  path  (source  to  interface  to  receiver)  is  not  very  close  to  the 
critical  angle.  If  this  path  is  close  to  critical,  then  the  reflection 
coefficient  magnitude  is  changing  rapidly  and  plane  wave  theory,  even  with 
displacement,  is  not  valid. 
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Whether  or  not  the  Stoneley  wave  is  favorably  excited  is  not  as  easy  to 
determine.  The  Stoneley  wave  pole  of  the  reflection  coefficient  must  be  found 
and  then  the  wave  number  corresponding  to  that  pole  must  be  used  in  equation  (31) 
to  establish  the  critical  geometry.  Since  we  only  need  the  geometry  approxi¬ 
mately,  it  should  be  acceptable  to  ignore  the  displacement  term  (the  right-hand- 
side)  of  equation  (31)  and,  also,  to  use  the  real  part  of  the  wave  number 
rather  than  the  true  saddle  point  value. 

To  conclude  this  section,  let  us  examine  several  examples  of  the  full 
wave  solution  and  compare  these  results  to  actual  measurements.  The  measure¬ 
ments,  made  by  Roever  and  Vining,10  are  high-quality  scale  model  measurements 
of  short  pulse  reflection  from  several  substances  with  both  compressional  and 
elastic  properties. 

Some  objections  may  be  raised  at  this  point  because  these  measurements 
are  scale-model  measurements  and  not  actual  ocean  measurements.  In  spite  of 
the  vast  quantity  of  ocean  acoustic  measurements  available  at  various  Navy 
and  university  laboratories,  measurements  that  are  well-defined  in  terms  of 
precise  experimental  geometry,  physical  parameters  of  the  medium,  and  charac¬ 
teristics  of  the  source  and  also  are  well  documented  as  to  experimental 
procedures  are  rare.  Initial  validation  of  a  model  against  the  usual  types 
of  measurements  would  involve  considerable  "adjustment"  of  the  physical 
parameters  in  order  to  get  good  agreement.  This  technique,  although  common, 
always  includes  the  possibility  that  the  medium  properties  were  adjusted  not  to 
correct  deficiencies  in  their  initial  estimates  but  to  offset  problems  in 
the  model  itself.  I  preferred,  in  this  investigation,  to  use  a  data  set  that 
I  could  model  without  such  manipulations. 

In  each  of  the  three  cases,  the  experiment  consisted  of  a  spark-gap  source 
and  a  barium  titanate  receiver  fixed  10  cm  apart  and  0.5  cm  above  the  reflecting 
material  in  salted  tap  water.  The  frequency  content  of  the  source  pulse  and 
the  medium  properties  and  experimental  geometry  were  such  that  the  horizontal 
range  was  between  10  and  0.1  A  (\  is  acoustic  wavelength).  This  presents  the 
most  difficult  modeling  problem  possible  as  wavefront  curvature  is  definitely 
important  and,  at  the  same  time,  the  asymptotic  methods  are  not  valid.  Hence, 
the  modeling  in  this  section  was  done  by  integration  of  the  field  integral, 
equation  (25),  and  subsequent  integration  over  frequency  by  FFT  using 
equation  (47). 

Physical  parameters  for  the  water  and  reflecting  media  were  measured  by 
Roever  and  Vining  and  included  compressional  and  shear  sound  speed,  compressional 
and  shear  attenuation,  and  density.  Since  these  properties  were  measured  on  the 
same  scale  as  the  experiment  itself,  the  property  estimates  were  used,  as 
given  in  reference  (10),  with  no  modifications. 

The  first  material  considered  was  a  high  viscosity  pitch  with  the  following 
acoustic  properties: 

Water  over  Pitch 

Cj  =  1496  m/sec,  c ^  =  2463  m/sec 

b  =  1003  m/sec,  .^A  j  =  1-27 
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A  =  1  dB/cm  at  500  kHz 
A$  =  0.6  dB/cm  at  25  kHz 

where 

=  sound  speed  in  water 

C2  =  congressional  speed  in  lower  medium 

b  =  shear  speed  in  lower  medium 

^  =  density  in  water 

$2  ~  density  in  lower  medium 

Ac  =  congressional  attenuation  in  lower  medium 

As  =  shear  attenuation  in  lower  medium. 

The  measured  and  the  predicted  time  waveforms  at  the  receiver  are  plotted  in 
figure  14. 

Overall,  the  agreement  between  the  predicted  signal  and  the  actual  measure 
ment  is  excellent.  The  first  feature,  arriving  at  about  50  vsec  after  trans¬ 
mission,  is  the  congressional  interface  wave.  Next,  just  prior  to  70  usee,  th; 
direct  and  simply-reflected  arrivals  appear  but  they  are  not  individually 
distinguishable  because  the  travel  time  difference  between  the  two  paths  is 
only  about  1/8  ^sec.  Finally,  the  Stoneley  wave  arrives  between  100  and 
150  usee.  Since  the  shear  sound  speed  in  the  pitch  is  less  than  the  sound 
speed  in  the  water,  there  is  no  critically  refracted  shear  wave  in  this  case. 

This  example  is  of  particular  practical  importance  because  ocean  sediments 
are  usually  only  weakly  elastic  with  low  shear  speed  and  high  shear  attenuation 
These  factors  might  lead  to  a  conclusion  that  the  shear  propagation  could  be 
neglected,  but  there  is  a  considerable  amount  of  energy  in  the  Stoneley  wave 
in  this  example.  As  long  as  the  reflector  has  some  shear  rigidity,  there  will 
be  a  Stoneley  wave  pole  in  the  reflection  coefficient.  In  this  case,  that 
pole  corresponds  to  a  phase  velocity  of  821  m/sec  which  is  considerably  lower 
than  the  shear  speed  of  the  pitch.  Hence,  the  Stoneley  wave  arrives  last. 

One  very  important  consideration  when  modeling  shear  propagation  in 
sediments  is  the  inclusion  of  shear  attenuation.  Since  the  Stoneley  wave  can 
be  a  fairly  efficient  transporter  of  energy  at  short  ranges,  if  attenuation 
is  neglected,  the  Stoneley  wave  can  appear  to  dominate  the  arrival.  The 
water-over-pi tch  example  was  run  with  zero  shear  attenuation  and,  while  the 
events  before  100  ..sec  remained  virtually  unchanged,  the  Stoneley  wave  arrival 
grew  both  in  its  positive  and  negative  cycles  to  the  point  where  the  negative 
peak  (at  about  130  ..sec)  was  more  than  three  times  the  amplitude  of  the  direct 
arrival  peak  at  68  ..sec.  Obviously,  this  would  have  resulted  in  a  gross 
over-estimate  of  the  role  of  shear  propagation.  As  it  is  in  figure  14,  the 
level  of  the  predicted  Stonely  wave  is  slightly  less  than  the  actual  event. 

This  difference  is  not  surprising  considering  the  difficulty  of  measuring  shear 
attenutation. 
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PITCH 


Figure  14.  Predicted  time  series  for  reflection  from  pitch  by  spherical  wave  theory. 
Measurement  is  displaced  below  predicted  curve. 
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The  next  reflecting  medium  was  cured  pi aster-of-pari s  and  the  experimentally 
determined  properties  were: 

Water  over  Plaster 

c^  =  1484  m/sec,  c2  =  3192  m/sec 

b  =  1814  m/sec,  e^0!  =  1-89 

Ac  =  0.6  dB/cm  at  500  kHz 

As  =  0.2  dB/cm  at  25  kHz. 

The  three  features  present  with  pitch  can  be  seen  here  again  in  figure  15. 

The  compressional  interface  wave  is  much  lower  in  level  (the  vertical  scale 
is  identical  to  that  of  figure  14)  appearing  at  40  ysec.  At  about  70  ysec, 
the  combined  direct/reflected  pulse  appears  and  a  rather  large  amplitude 
Stoneley  wave  appears  between  60  and  100  ysec.  The  direct  arrival  is  at  about 
the  same  place  as  before  because  the  water  has  almost  the  same  sound  speed  in 
both  cases,  but  both  the  compressional  interface  wave  and  the  Stoneley  wave 
have  arrived  earlier  because  the  speeds  in  plaster  are  higher  than  those  in 
pitch.  Also,  the  plaster  is  more  rigid  with  less  attenuation  of  shear  waves 
so  the  Stoneley  wave  amplitude  is  larger. 

An  additional  event  can  be  seen  at  60  ysec  and  this,  according  to  Strick,1® 
is  a  combination  of  a  small,  critically  refracted  shear  wave  and  a  psuedo- 
Rayleigh  wave.  Both  of  these  waves  are  interface  waves  that  travel  at  the 
shear  speed  of  the  lower  medium. 

Finally,  the  results  for  water  over  iron  are  shown  in  figure  16  and  the 
physical  properties  are  tabulated  below. 

Water  over  Iron 

c^  =  1500  m/sec,  =  5837  m/sec 

b  =  3247  m/sec,  ..^A  ^  =  7.87 

A  =  0.01  dB/cm  at  500  kHz 
c 

A$  =  0.01  dB/cm  at  25  kHz. 

No  attenuation  measurements  were  made  in  this  case  by  Roever  and  Vining  so 
very  small  values  were  used  merely  to  avoid  numerical  problems  associated 
with  integration  exactly  along  the  real  axis. 

Because  the  impedance  discontinuity  is  extremely  large  in  this  case  (the 
compressional  impedance  ratio  is  over  30),  the  combination  direct/reflected 
arrival  dominates  the  waveform.  The  Stoneley  wave  is  much  shorter  in  time 
than  before  and  only  the  negative  cycle  can  be  seen  at  about  66  „sec.  A  very 
small  compressional  interface  component  is  just  visible  at  25  -..sec  and,  again, 
there  is  a  combination  shear  interface  and  pseudo-Rayleigh  wave  at  40  usee. 

As  in  the  other  two  examples,  the  agreement  between  the  theory  and  the 
experiment  is  excellent. 
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waveform  is  lower  curve. 


IRON 


Figure  16.  Predicted  time  series  for  reflection  from  iron.  Lower  curve  is  measurement . 
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These  examples  are,  in  one  sense,  difficult  problems  because  the  wave 
effects  are  fully  developed  so  that  even  the  asymptotic  methods  do  not  work 
well.  They  are,  however,  simple  examples  in  that  the  reflector  is  a  homo¬ 
geneous  half-space.  More  complicated  reflector  structures  can  easily  be 
accommodated  though,  since  the  model  really  only  needs  a  reflection  coefficient 
and,  given  enough  time,  this  can  be  computed  for  very  complicated  structures 
(see  appendix  A). 


CONCLUSIONS 


In  this  work,  we  have  examined  two  general  types  of  models  for  propagation 
of  transient  acoustic  signals.  For  one  type,  we  developed  the  field  integral 
that  exactly  represents  the  acoustic  field  for  a  given  source  and  medium.  The 
other  type  is  actually  a  very  special  case  of  the  exact  solution.  If  enough 
assumptions  are  valid,  the  field  can  also  be  represented  by  ray  theory  with 
plane  wave  reflection  at  boundaries.  The  latter  model  is  much  simpler  and 
more  flexible;  however,  these  advantages  are  gained  by  sacrificing  some 
important  features  of  propagation  associated  with  wavefront  curvature.  In 
many  cases,  wavefront  curvature  is  negligible  but,  when  it  is  not,  some  form 
of  the  field  integral  must  be  evaluated. 

Whatever  model  is  used  must  suit  the  environment  in  question.  In  deep 
water  or  at  high  frequency,  the  plane  wave  model  may  be  quite  accurate.  At 
low  frequency,  near  the  ocean  floor  or  in  or  near  leaky  ducts,  the  plane  wave 

model  may  omit  significant  or  dominant  propagation  mechanisms. 

Fortunately,  there  are  some  techniques  that  are  intermediate  in  complexity 
and  still  offer  at  least  approximations  to  the  full  wave  solutions.  The 
saddle  point  or  steepest  descent  method  reduces  the  field  integral  to  one 
or  more  integrals  that  can  be  evaluated  over  short  intervals  and,  possibly, 
some  residues  or  branch  line  integrals.  As  long  as  the  true  saddle  points 
are  located  and  all  of  the  poles  and  branch  lines  are  accounted  for,  this 
approach  is  still  exact.  The  first  approximation  that  can  be  made  is  to 
integrate  only  over  a  small  region  near  each  saddle  point.  Also,  we  can 
neglect  the  offset  of  the  saddle  points  from  the  real  axis.  At  this  point, 
the  procedure  can  be  coupled  with  ray  theory  for  each  of  the  possible  paths 
and  some  representation  of  interface  waves,  trapped  modes,  and  field  dis¬ 
placement  will  be  maintained.  Finally,  if  the  phase  of  the  reflection 

coefficient  is  neglected  the  solution  reduces  to  classical  ray  theory  with 

a  single  specular  reflection  at  the  ocean  floor. 

The  decision  to  use  either  plane  wave  or  spherical  wave  methods  should 
be  based  on  the  wavefront  curvature  criterion.  If  an  acoustic  wavelength  is 
an  appreciable  part  of  the  ci rcumference  of  the  acoustic  wave  at  the  boundary, 
then  the  full  wave  solution  is  justified.  Some  weight  must  be  given  to  the 
relative  cost  of  the  approaches  though.  The  ray  theory  model  is  easily 
adapted  to  complex  multipath  situations  and  is  fast  and  inexpensive  to  run. 

In  some  cases,  this  may  justify  its  use  even  when  the  full  wave  solution  is 
indicated  by  the  wavefront  curvature  criterion,  as  long  as  the  shortcomings 
of  the  ray  approach  are  acknowledged. 
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We  have  seen  both  in  the  theoretical  development  and  in  the  experimental 
comparisons  that  the  models  can  be  used  to  interpret  propagation  results 
physically.  The  plane  wave  and  ray  theory  approach  obviously  treats  each 
path  distinctly.  The  saddle  point  method  also  breaks  the  solution  down  into 
a  number  of  physically  distinct  "paths"  and  special  wave  effects  may  be 
interpreted  in  light  of  these  components. 

Finally,  we  have  outlined  several  ways  in  which  the  FFT  can  be  used  to 
speed  the  calculations  required  by  these  models.  The  integration  over 
frequency  can  be  done  very  easily  using  the  FFT  and,  at  the  same  time,  the 
source  spectrum  can  be  included  so  that  the  resulting  response  is  not  the 
impulse  response  but,  instead,  a  simulation  of  the  complete  received  response. 
In  addition,  the  field  integration  over  wave  number  is  suited  to  evaluation 
by  FFT,  especially  if  results  are  required  at  many  different  source-to- 
receiver  ranges. 
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Appendix  A 

REFLECTION  FROM  MULTILAYERED  MEDIUM 


In  order  to  provide  for  future  expansion  of  the  reflection  model,  the  plane 
wave  reflection  coefficient  for  a  general,  multi-layered  medium  is  derived  below. 
Although  this  derivation  could  just  as  easily  be  done  for  any  quantity  that 
satisfies  the  wave  equation,  we  will  write  expressions  in  terms  of  velocity 
potential. 

The  generalized  medium  is  shown  in  figure  A-l.  In  each  layer  or  in  the 
fluid  half-space  above  the  layers,  the  expression  for  the  field  is, 

*  a1p1(z)  +  b1q1(z)  (zi_1  <  z  <  z^  (A-l) 


where 


a^ ,bj  =  constants 

p.j.q.j  =  independent  solutions  of  the  depth  equation 
«  velocity  potential  in  ith  layer 


The  equation  for  <f>0 ,  the  potential  in  the  fluid  half-space,  can  be  written  more 
specifically  in  terms  of  the  overall  reflection  coefficient  V^, 


=  const 


r-ig  z  is  z*l 

|f  0  °J 


-IB  z  i B_Z 

=  a  e  0  +  be 
o  o 


therefore. 


A-2) 


where 


ri  =  W 


The  boundary  conditions  between  each  pair  of  layers  are, 
(1)  continuity  of  pressure 

:i*i  (zi> 


A-l 
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(2)  continuity  of  normal  component  of  velocity 


ill 

8Z 


ali+l 


8Z 


zi 


2i 


These  two  equations  can  be  written  as  one  expression  relating  the  ratio  of 
coefficient  in  the  ith  layer  to  the  ratio  of  coefficients  in  the  layer  below. 


ri  = 


.  -^i  £Wi+i  +  +  fl'i  ^ri+iPi+i + 


(A-3) 


Pi  Cri+1pl+i  +  q!+1]  -  pi  [ri+1p1+1  +  qi+1]  Pi+1/"i 


where  the  p  and  q  function  are  evaluated  at  z= z-j  and  the  primes  indicate 
differentiation  with  respect  to  z. 

The  solution  process  generally  begins  at  the  solid  basement  half-space  and 
proceeds  upwards  layer  by  layer  until  r0  and,  therefore,  are  found.  Equation 
(A-3)  relates  any  two  adjacent  layers  and  equation  (A-2)  yields  the  total 
reflection  coefficient  Vt-  In  order  to  start  the  process  at  the  fluid-solid 
interface  at  zn,  we  use  the  impedance  condition,  equation  (5),  and  solve  for 
the  initial  coefficient  ratio  r^, 

r  .  -W  *  *  (V  (A.4) 

N  «n> 


where  the  bottom  impedance  involves  the  elastic  half-space  reflection  coefficient 
from  equation  (4) . 


If  the  sound  speed  in  the  layered  bottom  becomes  considerably  greater  than 
the  phase  velocity  of  the  incident  wavefront  (that  is,  the  vertex  velocity 
of  the  equivalent  ray),  then  it  is  neither  necessary  nor  desireable  to  start 
at  the  elastic  half-space.  As  the  medium  sound  speed  increases  beyond  the  phase 
velocity,  the  depth  function  decays  in  an  approximately  exponential  manner. 
Starting  the  solution  all  the  way  down  at  the  elastic  half-space  would  mean 
starting  with  an  extremely  small  field  amplitude.  Any  error  in  the  original 
value  would  grow  unacceptably  large  by  the  time  the  top  of  the  reflecting  layers 
was  reached. 


This  problem  is  avoided  Dy  including  a  false  basement  of  homogeneous  fluid 
at  some  depth  where  the  local  sound  speed  is  substantially  greater  than  the 
phase  velocity.  If  the  sound  speed  of  the  false  bottom  is  matched  to  the  real 
profile  at  this  point,  the  equivalent  impedance  is, 


i/.Ji 


A-3 
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where 


if  =  depth  of  false  bottom. 

This  impedance  value  is  inserted  into  equation  (A-4)  and  the  solution  progresses 
from  Zf  to  the  reflector  surface  at  z=0  in  the  usual  way. 


Appendix  B 
AIRY  FUNCTIONS 
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Appendix  B 
AIRY  FUNCTIONS 


Origin  of  Airy  Functions 

Airy  functions  arise  as  special  solutions  of  the  depth  equation  that  results 
from  separation  of  variables  of  the  wave  equation.  This  wave  equation,  in 
axial ly-symmetric  cylindrical  coordinates. 


I  JL_  (r^l\  +  III  -  1_  sff 

r  3r  \3r/  3z2  "  C2  3t2 

where 


c  =  sound  speed 
r  =  horizontal  range 
t  =  time 
z  =  depth 

v  =  field  quantity  (pressure,  for  example), 


reduces  to  a  Helmholtz  equation  for  harmonic  disturbances  and  then  separates 
into  a  range  equation  involving  Bessel  functions  and  a  depth  equation. 


+  [k2(z)  -  72]  U  =  0 

dz 


(B-l) 


where 


k  =  uj/c  ( z ) 

U  =  depth  function 

y  =  horizontal  wave  number  from  range  equation 
jj  =  angular  frequency  (2nf). 

If  the  sound  speed  c  is  constant,  the  solutions  of  equation  (B-l)  are 
sines  and  cosines;  however,  we  would  like  to  be  able  to  specify  some  variation 
in  sound  speed  with  depth.  Since  the  z-variation  of  equation  (B-l)  is  of  the 
form  l/c2,  we  can  make  the  equation  linear  in  z  if  we  let. 


~2 

c 


1 

(z) 


(B-2) 


where 

cQ  =  sound  speed  at  z  =  zQ 

g  - Tl/(z2  ’  zl^  =  pseudo-gradient. 

VI  c2  / 
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thus,  equation  (B-l)  can  be  written  as, 

^  +  ^  "  9Uj2(Z  "  Zq)  U  =  °'  (B"3) 
Being  linear  in  z,  this  equation  can  be  simplified  by  the  following  transformation. 


z  -  a  -  y2 j  -  gw2(z  -  zQ) 


wt.  _  i. 

dz  ~  'agu) 

Equation  (B-3)  then  becomes, 

a2(g»2)2  4  ♦  i  ZU  -  0 
dZ^  a 

which  is  Stoke's  equation. 


-  ZU  =  0 


(B-4) 


,  .  ,n  2 \ -2/ 3 
a  -  -(gw  ) 


Z  =  q(z  -  z  )  ~~2  ( ~~2  ~ 

q  \Co 


where 


q  =  (gw  ) 

Stoke's  equation  is  second  order;  therefore,  it  has  two  linearly  independent 
solutions.  A  particular  set  of  these  solutions  comprises  Airy  functions  or, 
in  symbolic  notation,  A i ( Z )  and  Bi(Z). 

Standard  Series  Forms  For  Airy  Functions^ 

The  standard  form  of  the  ascending  series  is  a  linear  combination  of  the 
ordinary  power  series  solutions  to  equation  (B-4). 


Ai(Z)  =  a 1 f ( Z )  -  a2 g(Z) 
Bi(Z)  =  *3  ra1f(Z)  +  a2g(Z) 


(B-6) 
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where 


al  =  Ai (0) 
a2  =  Ai'(O) 


f(2)  =  1  +  Z3  +  ^ 


76  1-4-7 

L  9! 


Z9  +  ... 


g(z)  =  z  +  Jr  z4  +  z7  + 


2-5-8 

10! 


(B-7) 


Z10  + 


These  series  are  valid  for  any  complex  Z  although  as  the  magnitude  of  Z  increases, 
the  number  of  terms  required  for  a  given  accuracy  also  increases.  For  this 
reason,  various  other  series  representations  have  been  derived. 

Asymptotic  series,  as  they  are  called,  are  series  that  are  valid  when  the 
argument  approaches  some  limiting  value  (infinity  in  this  case).  They  are  not 
necessarily  convergent  series  and  the  best  accuracy  is  obtained  by  summing 
until  the  terms  start  to  increase  in  magnitude.  For  the  Airy  functions,  these 
asymptotic  series  are  derived  by  applying  the  method  of  steepest  descent  to  the 
Airy  integrals!!.  The  derivation  procedure  is  unimportant  for  this  discussion 
so  we  will  merely  list  the  various  forms  and  their  ranges  of  applicability^: 

Definitions  - 


co  =  1 


216kk! 


do  =  1 


6k- 1  ck 


-  2  7 3/ 2 
"3 


(-ir  c. 


S.  =  ::  (-l)k  d  ;'k 

D  k=0  K 


se  =  ■  (-1: 

e  k=0 


::  (-1)K  d, 


::  (-1)K  c. 


-  2k- 1 


-2k- 1 


1 
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Exponential  forms  - 


Ai(Z) 

TT7  Z' 

^'Sa 

'arq  Z/  <  ~ 

A  i  ‘  ( Z ) 

1 

'  T T7 

zV;  sb 

/ arg  11  <  -t 

bi  ;z; 

v  \  V% 

'  Sa 

/arg  Z/  <  -i/ 3 

B  i 1  ( Z ) 

t  " 

'  Sb 

/arg  Z/  <  */3 

Qsci 1 latory 

forms  - 

redefine  ;  =  |-( 

-Z)3'2 

A  i  ( Z )  X  ( 

-Z)-'4  Ts 

i- 

sin;’  - 
e 

SQCOS;'j  /arg  Z/ 

>  2^/3 

A  •  ‘  '  Z ) 

1  (-Z)4 

Ts  cos 
l  ep 

- 1  +  S  sin;’|  /arg  Z/ 
op  J 

>  2-/3 

3’  'V  ~( 

[ 

S  cos-, 1 
e 

+  SQsin:'j  /arg  Z/ 

>  2-/3 

~  ■  •  ;  I 

(-z)u 

s  sin-  ’ 
.  ep 

-  Sopcos-'J  /arg  Z/ 

>  2-/3 

Special  oscillatory  forms  for  Bi  - 

,  2  /7  +  i  /3  3/2 

'"edeMne  =  ^  (Ze  ) 


(B- 10) 


(B-ll) 


(B-12) 


(B-13) 
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Bi(Z)  =  e±i7r/4  Z_i 


-Socos(?'  +  i  ln2)]  /(arg  Z)+  j  /  <  2tt/3 

Bi'(Z)  -V7e+iTT/4  z*  [sepcos  (c‘  TJln2)  (B  1 

+  SQpsi n( c '  +  j  1 n2 )J  /(arg  Z)+  J  /  <  2i/3 

As  long  as  the  magnitude  of  Z  is  great  enough  to  insure  sufficient  accuracy  of 
the  asymptotic  series,  one  or  another  of  these  can  be  used  to  compute  the  Airy 
functions  or  their  derivatives  in  any  region  of  the  complex  Z-plane. 


Problems  In  Numerical  Computation 

It  is  evident  from  equations  (B-6)  and  (B-7)  that  there  is  a  region  on  and 
around  the  positive  real  Z-axis  in  which  the  terms  of  the  series  for  Ai  (and  Ai ' ) 
alternate  in  sign.  The  individual  terms  become  very  large  as  the  series 
progresses  thus  requiring  the  differences  of  large  numbers.  Since  these 
differences  should  be  small  as  the  series  converges,  the  numerical  range  imposed 
by  quantization  in  the  computer  word  is  exceeded  very  quickly. 

Consequently,  the  ascending  series  cannot  be  computed  far  enough  in  Z  to 
match  satisfactorily  with  the  asymptotic  series  of  equations  (B-10).  Even  with 
the  120-bit  double  precision  words  of  the  CDC  Cyber  170  computer,  only  seven 
decimal  digits  are  correct  at  the  best  matching  point  between  the  ascending 
and  asymptotic  series.  The  more  common  64-bit  double  precision  word  would 
restrict  the  accuracy  even  further,  perhaps  to  the  point  of  creating  problems 
in  iterative  calculations. 


Origin-Translated  Ascending  Series 

The  solution  to  this  problem  lies  in  modifying  the  ascending  series  so  that 
the  terms  do  not  alternate  in  sign.  This  is  easily  done  by  solving  the 
differential  equation  (B-4)  with  a  power  series  expanded  around  a  point  dis¬ 
placed  along  the  positive  real  Z-axis.  Let, 

t  =  Z  -  Z  (B-15) 

o 


where 

lQ  =  origin  point  for  the  series, 
so  that  equation  (B-l)  becomes, 

^7  -  (t  +  ZJ  U  =  0. 

dr  0 


(B-16) 


B-5 
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Set  up  the  power  series  solution  in  terms  of  the  new  variable  t, 


U  =  z  antn  ( B- 1 7 ) 

n=Q  n 


and  substitute  into  equation  (B-16).  By  equating  the  coefficients  of  like 
powers  of  t,  we  can  form  the  following  relationships  in  which  a  and  a,  are 
arbitrary,  0  1 


.  Vl  *  Zoa„ 

n+2  ( n+2 )  (n+l) 


(B-18) 


or 


U(Z)  =  +  axt  + 


V*  +  Ji  [ao  +  Zoal 


1  t3  + 


from  which  we  see  that, 
ao  ■  U<2o> 

(B-19) 


The  conditions  of  equations  (B-19)  allow  us  to  develop  the  Airy  functions  if  we 
know  Ai(Z0)  and  Ai’(Z0)  (or,  of  course,  Bi(Z0)  and  Bi'(Z0)).  Notice  that  Z0 
is  not  restricted  to  being  real  although  we  will  only  consider  the  real  case 
since  this  will  solve  our  particular  accuracy  problem. 

In  the  following  section,  we  will  consider  the  problem  of  finding  the 
starting  values  Ai(Z0)  and  Ai'(Zp)  as  well  as  the  errors  of  each  of  the  other 
series  over  their  regions  of  validity  in  the  Z-plane.  At  this  point,  it  should 
be  mentioned  that  the  origin-shifted  series  completely  eliminates  the  computa¬ 
tional  problems  involved  with  the  usual  ascending  series.  When  the  origin- 
shifted  series  is  evaluated  for  Z  less  than  Z0,  the  terms  are  no  longer 
alternating  in  sign  and  the  series  converges  rapidly. 

Error  Analysis 

Except  for  the  problem  region  (roughly  a  +30°  wedge  centered  on  the  positive 
real  Z-axis)  discussed  above,  the  standard  against  which  each  of  the  series  was 
evaluated  is  a  double  precision  (120-bit  word  or  about  29  decimal  digits) 
ascending  series  based  on  equations  (B-6).  After  some  experimentation,  I 
selected  /!/  =  6  as  the  crossover  point  between  the  single  precision  ascending 
series  and  the  asymptotic  series.  Figure  B-l  shows  the  number  of  digits  that 
correspond  exactly  with  the  double  precision  ascending  series  for  each  form  of 
the  single  precision  series  at  /!/  -  6.  The  problem  around  the  positive  real 
axis  can  be  seen  readily  in  the  exponential  case  and  the  ascending  series  case 


B-6 
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for  Ai  (note  the  circled  numbers).  In  addition,  we  can  see  over  what  regions 
the  other  series  should  be  used  if,  as  in  this  case,  we  would  like  at  least 
eight  decimal  digits  to  be  correct  using  only  single  precision  (60-bit). 

Since  no  series  shown  in  figure  B-l  covers  the  /arg  Z/  <  u/6  region 
adequately,  we  must  now  consider  the  origin-shifted  series.  The  error  of  an 
asymptotic  series  is  shown  in  reference  (11)  to  be  of  the  order  of  the  last 

computed  because  beyond  this  point  the 
successive  termsl2, 


the  it*1  term  or. 


c;  therefore,  for  any  g,  this  point  can  be 
ith  term  then  gives  the  order  of  accuracy 

of  the  summation. 

In  this  work,  the  exponential  asymptotic  series  gave  a  satisfactory  v  lue 
at  Z  =  12,  so  I  first  generated  an  origin-shifted  series  about  this  point.  As 
shown  in  figure  2,  this  series  tracked  well  with  the  asymptotic  series  down 
to  Z  =  7.  Next,  I  constructed  a  series  with  its  origin  at  Z  =  8  using  the 
values  of  Ai(8)  and  Ai ' (8)  from  the  Z  =  12  series.  This  new  series  was  evalu¬ 
ated  all  the  way  back  to  Z  =  0  and,  as  seen  in  figure  B-2,  the  agreement  with 
the  double  precision  zero-origin  series  is  excellent  below  Z  =  4.  This  shows 
that  even  the  double  precision  (120-bit)  ascending  series  deteriorates  rapidly 
beyond  Z  =  5.  Finally,  the  Z  =  8  series  was  used  to  compute  Ai ( 6 )  and  Ai * (6 ) 
the  starter  values  for  the  series  actually  used  in  the  Airy  function  algorithm 
These  values  could  have  been  calculated  from  the  Z  =  12  series  to  roughly  the 
same  accuracy  but  this  was  not  apparent  until  the  Z  =  8  series  had  been  tested 

Routine  For  Calculation  of  Airy  Functions 

Based  on  the  error  study  summarized  in  the  last  section,  the  argument 
plane  (the  complex  Z-plane)  was  divided  into  the  regions  shown  in  figure  B-3. 
These  regions  and  the  associated  forms  of  the  Airy  function  series  are  as 
f o 1 1 ows , 

I)  HI  <  3  :  ascending  series  (Z0  =  0)  for  Ai  and  Bi  -  equations  (B-6), 

II)  3  <_ /Z  /  <_6:  ascending  series  (Z0  =  0)  for  Bi  -  equation  (B-6)  and 

a)  /arg  Z/  ^  ~/6;  ascending  series  (Z  =  6)  for  Ai  - 
equation  (B-17), 

b)  /arg  Z/  >  ;/6:  ascending  series  (Z0  =  0)  for  Ai  - 
equation  (B-6), 


term  included.  This  term  is  easily 
terms  grow.  Hence,  we  must  compare 


V 


-1 


ai+l? 


-i-1 


until  the  i+lth  term  is  larger  than 


Vi 


>  1 


The  coefficients  are  independent  of 
easily  calculated.  The  size  of  the 
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Ill)  /Z/  >_6:  a)  /argZ/  >_  2tt/3:  oscillatory  asymptotic  series  for  Ai 
and  Bi  -  equations  (B-12)  and  (B-13), 

b)  /arg  Z/  <  2n/3:  exponential  series  for  Ai  -  equation 
(B-10)  and  special  oscillatory  series  for  Bi  - 
equation  (B-14). 

The  starting  values  for  the  ascending  series  are: 

I) 6  z0  *  0 

Ai (0)  »  0.355  028  053  888 

Ai'(0)  =  -0.258  819  403  793 

Bi (0)  =  Ai  (0) 

Bi 1 (0)  =  -/3  Ai'(0) 

II)  Zq=6 

Ai (6)  =  0.994  769  436  025  x  10‘5 

Ai 1 (6)  =  -0.247  652  003  970  x  10“4 


8-8 


Number  of  correct  digits  for  single  precision  (60-bit)  series.  Standard  for  comparison  is  double  precision 
ascending  series.  Circles  numbers  indicate  breakdown  of  this  standard  and  are  not  otherwise  meaningful. 
Outer  ring  refers  to  Bi  and  Bi '  for  /Z /  =  6  and  0  <_  argZ  tt  in  tt/6  steps.  Inner  ring  refers  to  Ai  and  Ai ! 
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